Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 <cmath>
61#include <TMath.h>
62#include <TRandom3.h>
63#include <TFile.h>
64#include <TTree.h>
65#include <TLorentzVector.h>
66
67#define MASS1 0.511E-3
68
69using std::cout;
70
71TRandom *g_rnd=nullptr;
72
73class ToyEvent7 {
74public:
75 void GenerateDataEvent(TRandom *rnd);
76 void GenerateSignalEvent(TRandom *rnd);
77 void GenerateBgrEvent(TRandom *rnd);
78 // reconstructed quantities
79 inline Double_t GetMRec(int i) const { return fMRec[i]; }
80 inline Double_t GetPtRec(int i) const { return fPtRec[i]; }
81 inline Double_t GetEtaRec(int i) const { return fEtaRec[i]; }
82 inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
83 inline Double_t GetPhiRec(int i) const { return fPhiRec[i]; }
84 inline Bool_t IsTriggered(void) const { return fIsTriggered; }
85
86 // generator level quantities
87 inline Double_t GetMGen(int i) const {
88 if(IsSignal()) return fMGen[i];
89 else return -1.0;
90 }
91 inline Double_t GetPtGen(int i) const {
92 if(IsSignal()) return fPtGen[i];
93 else return -1.0;
94 }
95 inline Double_t GetEtaGen(int i) const {
96 if(IsSignal()) return fEtaGen[i];
97 else return 999.0;
98 }
99 inline Double_t GetPhiGen(int i) const {
100 if(IsSignal()) return fPhiGen[i];
101 else return 999.0;
102 }
103 inline Bool_t IsSignal(void) const { return fIsSignal; }
104protected:
105
106 void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
107 void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
108 void GenerateReco(TRandom *rnd);
109
110 // reconstructed quantities
111 Double_t fMRec[3];
112 Double_t fPtRec[3];
113 Double_t fEtaRec[3];
114 Double_t fPhiRec[3];
115 Double_t fDiscriminator;
116 Bool_t fIsTriggered;
117 // generated quantities
118 Double_t fMGen[3];
119 Double_t fPtGen[3];
120 Double_t fEtaGen[3];
121 Double_t fPhiGen[3];
122 Bool_t fIsSignal;
123public:
124 static Double_t kDataSignalFraction;
125 static Double_t kMCSignalFraction;
126
127};
128
129void testUnfold7a()
130{
131 // random generator
132 g_rnd=new TRandom3(4711);
133
134 // data and MC number of events
135 Double_t muData0=5000.;
136 // luminosity error
137 Double_t muData=muData0*g_rnd->Gaus(1.0,0.03);
138 // stat error
139 Int_t neventData = g_rnd->Poisson( muData);
140
141 // generated number of MC events
142 Int_t neventSigmc = 250000;
143 Int_t neventBgrmc = 100000;
144
145 Float_t etaRec[3],ptRec[3],phiRec[3],mRec[3],discr;
146 Float_t etaGen[3],ptGen[3],phiGen[3],mGen[3];
147 Float_t weight;
148 Int_t istriggered,issignal;
149
150 //==================================================================
151 // Step 1: generate data TTree
152
153 TFile *dataFile=new TFile("testUnfold7_data.root","recreate");
154 TTree *dataTree=new TTree("data","event");
155
156 dataTree->Branch("etarec",etaRec,"etarec[3]/F");
157 dataTree->Branch("ptrec",ptRec,"ptrec[3]/F");
158 dataTree->Branch("phirec",phiRec,"phirec[3]/F");
159 dataTree->Branch("mrec",mRec,"mrec[3]/F");
160 dataTree->Branch("discr",&discr,"discr/F");
161
162 // for real data, only the triggered events are available
163 dataTree->Branch("istriggered",&istriggered,"istriggered/I");
164 // data truth parameters
165 dataTree->Branch("etagen",etaGen,"etagen[3]/F");
166 dataTree->Branch("ptgen",ptGen,"ptgen[3]/F");
167 dataTree->Branch("phigen",phiGen,"phigen[3]/F");
168 dataTree->Branch("mgen",mGen,"mgen[3]/F");
169 dataTree->Branch("issignal",&issignal,"issignal/I");
170
171 cout<<"fill data tree\n";
172
173 //Int_t nEvent=0,nTriggered=0;
174 for(int ievent=0;ievent<neventData;ievent++) {
175 ToyEvent7 event;
176 event.GenerateDataEvent(g_rnd);
177 for(int i=0;i<3;i++) {
178 etaRec[i]=event.GetEtaRec(i);
179 ptRec[i]=event.GetPtRec(i);
180 phiRec[i]=event.GetPhiRec(i);
181 mRec[i]=event.GetMRec(i);
182 etaGen[i]=event.GetEtaGen(i);
183 ptGen[i]=event.GetPtGen(i);
184 phiGen[i]=event.GetPhiGen(i);
185 mGen[i]=event.GetMGen(i);
186 }
187 discr=event.GetDiscriminator();
188 istriggered=event.IsTriggered() ? 1 : 0;
189 issignal=event.IsSignal() ? 1 : 0;
190
191 dataTree->Fill();
192
193 if(!(ievent%100000)) cout<<" data event "<<ievent<<"\n";
194
195 //if(istriggered) nTriggered++;
196 //nEvent++;
197
198 }
199
200 dataTree->Write();
201 delete dataTree;
202 delete dataFile;
203
204 //==================================================================
205 // Step 2: generate signal TTree
206
207 TFile *signalFile=new TFile("testUnfold7_signal.root","recreate");
208 TTree *signalTree=new TTree("signal","event");
209
210 signalTree->Branch("etarec",etaRec,"etarec[3]/F");
211 signalTree->Branch("ptrec",ptRec,"ptrec[3]/F");
212 signalTree->Branch("phirec",ptRec,"phirec[3]/F");
213 signalTree->Branch("mrec",mRec,"mrec[3]/F");
214 signalTree->Branch("discr",&discr,"discr/F");
215
216 // for real data, only the triggered events are available
217 signalTree->Branch("istriggered",&istriggered,"istriggered/I");
218 // data truth parameters
219 signalTree->Branch("etagen",etaGen,"etagen[3]/F");
220 signalTree->Branch("ptgen",ptGen,"ptgen[3]/F");
221 signalTree->Branch("phigen",phiGen,"phigen[3]/F");
222 signalTree->Branch("weight",&weight,"weight/F");
223 signalTree->Branch("mgen",mGen,"mgen[3]/F");
224
225 cout<<"fill signal tree\n";
226
227 weight=ToyEvent7::kMCSignalFraction*muData0/neventSigmc;
228
229 for(int ievent=0;ievent<neventSigmc;ievent++) {
230 ToyEvent7 event;
231 event.GenerateSignalEvent(g_rnd);
232
233 for(int i=0;i<3;i++) {
234 etaRec[i]=event.GetEtaRec(i);
235 ptRec[i]=event.GetPtRec(i);
236 phiRec[i]=event.GetPhiRec(i);
237 mRec[i]=event.GetMRec(i);
238 etaGen[i]=event.GetEtaGen(i);
239 ptGen[i]=event.GetPtGen(i);
240 phiGen[i]=event.GetPhiGen(i);
241 mGen[i]=event.GetMGen(i);
242 }
243 discr=event.GetDiscriminator();
244 istriggered=event.IsTriggered() ? 1 : 0;
245
246 if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
247
248 signalTree->Fill();
249 }
250
251 signalTree->Write();
252 delete signalTree;
253 delete signalFile;
254
255 // ==============================================================
256 // Step 3: generate background MC TTree
257
258 TFile *bgrFile=new TFile("testUnfold7_background.root","recreate");
259 TTree *bgrTree=new TTree("background","event");
260
261 bgrTree->Branch("etarec",&etaRec,"etarec[3]/F");
262 bgrTree->Branch("ptrec",&ptRec,"ptrec[3]/F");
263 bgrTree->Branch("phirec",&phiRec,"phirec[3]/F");
264 bgrTree->Branch("mrec",&mRec,"mrec[3]/F");
265 bgrTree->Branch("discr",&discr,"discr/F");
266 bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
267 bgrTree->Branch("weight",&weight,"weight/F");
268
269 cout<<"fill background tree\n";
270
271 weight=(1.-ToyEvent7::kMCSignalFraction)*muData0/neventBgrmc;
272
273 for(int ievent=0;ievent<neventBgrmc;ievent++) {
274 ToyEvent7 event;
275 event.GenerateBgrEvent(g_rnd);
276 for(int i=0;i<3;i++) {
277 etaRec[i]=event.GetEtaRec(i);
278 ptRec[i]=event.GetPtRec(i);
279 phiRec[i]=event.GetPhiRec(i);
280 }
281 discr=event.GetDiscriminator();
282 istriggered=event.IsTriggered() ? 1 : 0;
283
284 if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
285
286 bgrTree->Fill();
287 }
288
289 bgrTree->Write();
290 delete bgrTree;
291 delete bgrFile;
292}
293
294Double_t ToyEvent7::kDataSignalFraction=0.75;
295Double_t ToyEvent7::kMCSignalFraction=0.75;
296
297void ToyEvent7::GenerateDataEvent(TRandom *rnd) {
298 fIsSignal=rnd->Uniform()<kDataSignalFraction;
299 if(IsSignal()) {
300 GenerateSignalKinematics(rnd,kTRUE);
301 } else {
302 GenerateBgrKinematics(rnd,kTRUE);
303 }
304 GenerateReco(rnd);
305}
306
307void ToyEvent7::GenerateSignalEvent(TRandom *rnd) {
308 fIsSignal=true;
309 GenerateSignalKinematics(rnd,kFALSE);
310 GenerateReco(rnd);
311}
312
313void ToyEvent7::GenerateBgrEvent(TRandom *rnd) {
314 fIsSignal=false;
315 GenerateBgrKinematics(rnd,kFALSE);
316 GenerateReco(rnd);
317}
318
319void ToyEvent7::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
320
321 // fake decay of Z0 to two fermions
322 double M0=91.1876;
323 double Gamma=2.4952;
324 // generated mass
325 do {
326 fMGen[2]=rnd->BreitWigner(M0,Gamma);
327 } while(fMGen[2]<=0.0);
328
329 double N_ETA=3.0;
330 double MU_PT=5.;
331 double SIGMA_PT=2.0;
332 double DECAY_A=0.2;
333 if(isData) {
334 //N_ETA=2.5;
335 MU_PT=6.;
336 SIGMA_PT=1.8;
337 //DECAY_A=0.5;
338 }
339 fEtaGen[2]=TMath::Power(rnd->Uniform(0,1.5),N_ETA);
340 if(rnd->Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
341 fPhiGen[2]=rnd->Uniform(-M_PI,M_PI);
342 do {
343 fPtGen[2]=rnd->Landau(MU_PT,SIGMA_PT);
344 } while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
345 //========================== decay
347 sum.SetPtEtaPhiM(fPtGen[2],fEtaGen[2],fPhiGen[2],fMGen[2]);
348 // boost into lab-frame
349 TVector3 boost=sum.BoostVector();
350 // decay in rest-frame
351
352 TLorentzVector p[3];
353 double m=MASS1;
354 double costh;
355 do {
356 double r=rnd->Uniform(-1.,1.);
357 costh=r*(1.+DECAY_A*r*r);
358 } while(fabs(costh)>=1.0);
359 double phi=rnd->Uniform(-M_PI,M_PI);
360 double e=0.5*sum.M();
361 double ptot=TMath::Sqrt(e+m)*TMath::Sqrt(e-m);
362 double pz=ptot*costh;
363 double pt=TMath::Sqrt(ptot+pz)*TMath::Sqrt(ptot-pz);
364 double px=pt*cos(phi);
365 double py=pt*sin(phi);
366 p[0].SetXYZT(px,py,pz,e);
367 p[1].SetXYZT(-px,-py,-pz,e);
368 for(int i=0;i<2;i++) {
369 p[i].Boost(boost);
370 }
371 p[2]=p[0]+p[1];
372 for(int i=0;i<3;i++) {
373 fPtGen[i]=p[i].Pt();
374 fEtaGen[i]=p[i].Eta();
375 fPhiGen[i]=p[i].Phi();
376 fMGen[i]=p[i].M();
377 }
378}
379
380void ToyEvent7::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
381 for(int i=0;i<3;i++) {
382 fPtGen[i]=0.0;
383 fEtaGen[i]=0.0;
384 fPhiGen[i]=0.0;
385 }
386 TLorentzVector p[3];
387 for(int i=0;i<2;i++) {
388 p[i].SetPtEtaPhiM(rnd->Exp(15.0),rnd->Uniform(-3.,3.),
389 rnd->Uniform(-M_PI,M_PI),isData ? MASS1 : MASS1);
390 }
391 p[2]=p[0]+p[1];
392 for(int i=0;i<3;i++) {
393 fPtRec[i]=p[i].Pt();
394 fEtaRec[i]=p[i].Eta();
395 fPhiRec[i]=p[i].Phi();
396 fMRec[i]=p[i].M();
397 }
398}
399
400void ToyEvent7::GenerateReco(TRandom *rnd) {
401 if(fIsSignal) {
402 TLorentzVector p[3];
403 for(int i=0;i<2;i++) {
404 Double_t expEta=TMath::Exp(fEtaGen[i]);
405 Double_t coshEta=(expEta+1./expEta);
406 Double_t eGen=fPtGen[i]*coshEta;
407 Double_t sigmaE=
408 0.1*TMath::Sqrt(eGen)
409 +1.0*coshEta
410 +0.01*eGen;
411 Double_t eRec;
412 do {
413 eRec=rnd->Gaus(eGen,sigmaE);
414 } while(eRec<=0.0);
415 Double_t sigmaEta=0.1+0.02*TMath::Abs(fEtaGen[i]);
416 p[i].SetPtEtaPhiM(eRec/(expEta+1./expEta),
417 rnd->Gaus(fEtaGen[i],sigmaEta),
418 remainder(rnd->Gaus(fPhiGen[i],0.03),2.*M_PI),
419 MASS1);
420 }
421 p[2]=p[0]+p[1];
422 for(int i=0;i<3;i++) {
423 fPtRec[i]=p[i].Pt();
424 fEtaRec[i]=p[i].Eta();
425 fPhiRec[i]=p[i].Phi();
426 fMRec[i]=p[i].M();
427 }
428 }
429 if(fIsSignal) {
430 do {
431 Double_t tauDiscr=0.08-0.04/(1.+fPtRec[2]/10.0);
432 Double_t sigmaDiscr=0.01;
433 fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
434 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
435 } else {
436 do {
437 Double_t tauDiscr=0.15-0.05/(1.+fPtRec[2]/5.0)+0.1*fEtaRec[2];
438 Double_t sigmaDiscr=0.02+0.01*fEtaRec[2];
439 fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
440 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
441 }
442 fIsTriggered=false;
443 for(int i=0;i<2;i++) {
444 if(rnd->Uniform()<0.92/(TMath::Exp(-(fPtRec[i]-15.5)/2.5)+1.)) fIsTriggered=true;
445 }
446}
#define e(i)
Definition RSha256.hxx:103
#define M_PI
Definition Rotated.cxx:105
bool Bool_t
Definition RtypesCore.h:63
int Int_t
Definition RtypesCore.h:45
float Float_t
Definition RtypesCore.h:57
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
Definition TFile.h:53
void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m)
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
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:404
virtual Double_t Exp(Double_t tau)
Returns an exponential deviate.
Definition TRandom.cxx:252
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:682
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:381
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:226
A TTree represents a columnar dataset.
Definition TTree.h:79
virtual Int_t Fill()
Fill all branches.
Definition TTree.cxx:4603
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:353
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) override
Write this object to the current directory.
Definition TTree.cxx:9743
TPaveText * pt
RVec< PromoteTypes< T0, T1 > > remainder(const T0 &x, const RVec< T1 > &v)
Definition RVec.hxx:1798
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition RVec.hxx:1815
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition RVec.hxx:1814
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:366
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Definition TMath.h:709
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345