Logo ROOT   6.08/07
Reference Guide
testUnfold5a.C
Go to the documentation of this file.
1 /// \defgroup tutorial_unfold5 TUnfoldDensity and TUnfoldBinning test suite
2 /// \ingroup tutorial_unfold
3 ///
4 /// This is an example of unfolding a two-dimensional distribution.
5 /// Also using an auxiliary measurement to constrain some background
6 ///
7 /// The example comprises several macros:
8 ///
9 /// - testUnfold5a.C create root files with TTree objects for
10 /// signal, background and data.
11 /// - write files:
12 /// - testUnfold5_signal.root
13 /// - testUnfold5_background.root
14 /// - testUnfold5_data.root
15 ///
16 /// - testUnfold5b.C create a root file with the TUnfoldBinning objects
17 /// -> write file testUnfold5_binning.root
18 ///
19 /// - testUnfold5c.C loop over trees and fill histograms based on the
20 /// TUnfoldBinning objects
21 /// - read:
22 /// - testUnfold5_binning.root
23 /// - testUnfold5_signal.root
24 /// - testUnfold5_background.root
25 /// - testUnfold5_data.root
26 /// - write:
27 /// - testUnfold5_histograms.root
28 ///
29 /// - testUnfold5d.C run the unfolding
30 /// - read:
31 /// - testUnfold5_histograms.root
32 /// - write:
33 /// - testUnfold5_result.root
34 /// - testUnfold5_result.ps
35 ///
36 /// \file
37 /// \ingroup tutorial_unfold5
38 /// \notebook -nodraw
39 /// Version 17.0 example for multi-dimensional unfolding
40 ///
41 /// \macro_output
42 /// \macro_code
43 ///
44 /// \author Stefan Schmitt, DESY
45 
46 #include <iostream>
47 #include <map>
48 #include <cmath>
49 #include <TMath.h>
50 #include <TRandom3.h>
51 #include <TFile.h>
52 #include <TTree.h>
53 
54 using namespace std;
55 
56 TRandom *g_rnd=0;
57 
58 class ToyEvent {
59 public:
60  void GenerateDataEvent(TRandom *rnd);
61  void GenerateSignalEvent(TRandom *rnd);
62  void GenerateBgrEvent(TRandom *rnd);
63  // reconstructed quantities
64  inline Double_t GetPtRec(void) const { return fPtRec; }
65  inline Double_t GetEtaRec(void) const { return fEtaRec; }
66  inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
67  inline Bool_t IsTriggered(void) const { return fIsTriggered; }
68 
69  // generator level quantities
70  inline Double_t GetPtGen(void) const {
71  if(IsSignal()) return fPtGen;
72  else return -1.0;
73  }
74  inline Double_t GetEtaGen(void) const {
75  if(IsSignal()) return fEtaGen;
76  else return 999.0;
77  }
78  inline Bool_t IsSignal(void) const { return fIsSignal; }
79 protected:
80 
81  void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
82  void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
83  void GenerateReco(TRandom *rnd);
84 
85  // reconstructed quantities
86  Double_t fPtRec;
87  Double_t fEtaRec;
88  Double_t fDiscriminator;
89  Bool_t fIsTriggered;
90  // generated quantities
91  Double_t fPtGen;
92  Double_t fEtaGen;
93  Bool_t fIsSignal;
94 
95  static Double_t kDataSignalFraction;
96 
97 };
98 
99 void testUnfold5a()
100 {
101  // random generator
102  g_rnd=new TRandom3();
103 
104  // data and MC number of events
105  const Int_t neventData = 20000;
106  Double_t const neventSignalMC =2000000;
107  Double_t const neventBgrMC =2000000;
108 
109  Float_t etaRec,ptRec,discr,etaGen,ptGen;
110  Int_t istriggered,issignal;
111 
112  //==================================================================
113  // Step 1: generate data TTree
114 
115  TFile *dataFile=new TFile("testUnfold5_data.root","recreate");
116  TTree *dataTree=new TTree("data","event");
117 
118  dataTree->Branch("etarec",&etaRec,"etarec/F");
119  dataTree->Branch("ptrec",&ptRec,"ptrec/F");
120  dataTree->Branch("discr",&discr,"discr/F");
121 
122  // for real data, only the triggered events are available
123  dataTree->Branch("istriggered",&istriggered,"istriggered/I");
124  // data truth parameters
125  dataTree->Branch("etagen",&etaGen,"etagen/F");
126  dataTree->Branch("ptgen",&ptGen,"ptgen/F");
127  dataTree->Branch("issignal",&issignal,"issignal/I");
128 
129  cout<<"fill data tree\n";
130 
131  Int_t nEvent=0,nTriggered=0;
132  while(nTriggered<neventData) {
133  ToyEvent event;
134  event.GenerateDataEvent(g_rnd);
135 
136  etaRec=event.GetEtaRec();
137  ptRec=event.GetPtRec();
138  discr=event.GetDiscriminator();
139  istriggered=event.IsTriggered() ? 1 : 0;
140  etaGen=event.GetEtaGen();
141  ptGen=event.GetPtGen();
142  issignal=event.IsSignal() ? 1 : 0;
143 
144  dataTree->Fill();
145 
146  if(!(nEvent%100000)) cout<<" data event "<<nEvent<<"\n";
147 
148  if(istriggered) nTriggered++;
149  nEvent++;
150 
151  }
152 
153  dataTree->Write();
154  delete dataTree;
155  delete dataFile;
156 
157  //==================================================================
158  // Step 2: generate signal TTree
159 
160  TFile *signalFile=new TFile("testUnfold5_signal.root","recreate");
161  TTree *signalTree=new TTree("signal","event");
162 
163  signalTree->Branch("etarec",&etaRec,"etarec/F");
164  signalTree->Branch("ptrec",&ptRec,"ptrec/F");
165  signalTree->Branch("discr",&discr,"discr/F");
166  signalTree->Branch("istriggered",&istriggered,"istriggered/I");
167  signalTree->Branch("etagen",&etaGen,"etagen/F");
168  signalTree->Branch("ptgen",&ptGen,"ptgen/F");
169 
170  cout<<"fill signal tree\n";
171 
172  for(int ievent=0;ievent<neventSignalMC;ievent++) {
173  ToyEvent event;
174  event.GenerateSignalEvent(g_rnd);
175 
176  etaRec=event.GetEtaRec();
177  ptRec=event.GetPtRec();
178  discr=event.GetDiscriminator();
179  istriggered=event.IsTriggered() ? 1 : 0;
180  etaGen=event.GetEtaGen();
181  ptGen=event.GetPtGen();
182 
183  if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
184 
185  signalTree->Fill();
186  }
187 
188  signalTree->Write();
189  delete signalTree;
190  delete signalFile;
191 
192  // ==============================================================
193  // Step 3: generate background MC TTree
194 
195  TFile *bgrFile=new TFile("testUnfold5_background.root","recreate");
196  TTree *bgrTree=new TTree("background","event");
197 
198  bgrTree->Branch("etarec",&etaRec,"etarec/F");
199  bgrTree->Branch("ptrec",&ptRec,"ptrec/F");
200  bgrTree->Branch("discr",&discr,"discr/F");
201  bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
202 
203  cout<<"fill background tree\n";
204 
205  for(int ievent=0;ievent<neventBgrMC;ievent++) {
206  ToyEvent event;
207  event.GenerateBgrEvent(g_rnd);
208  etaRec=event.GetEtaRec();
209  ptRec=event.GetPtRec();
210  discr=event.GetDiscriminator();
211  istriggered=event.IsTriggered() ? 1 : 0;
212 
213  if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
214 
215  bgrTree->Fill();
216  }
217 
218  bgrTree->Write();
219  delete bgrTree;
220  delete bgrFile;
221 }
222 
223 Double_t ToyEvent::kDataSignalFraction=0.8;
224 
225 void ToyEvent::GenerateDataEvent(TRandom *rnd) {
226  fIsSignal=rnd->Uniform()<kDataSignalFraction;
227  if(IsSignal()) {
228  GenerateSignalKinematics(rnd,kTRUE);
229  } else {
230  GenerateBgrKinematics(rnd,kTRUE);
231  }
232  GenerateReco(rnd);
233 }
234 
235 void ToyEvent::GenerateSignalEvent(TRandom *rnd) {
236  fIsSignal=1;
237  GenerateSignalKinematics(rnd,kFALSE);
238  GenerateReco(rnd);
239 }
240 
241 void ToyEvent::GenerateBgrEvent(TRandom *rnd) {
242  fIsSignal=0;
243  GenerateBgrKinematics(rnd,kFALSE);
244  GenerateReco(rnd);
245 }
246 
247 void ToyEvent::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
248  Double_t e_T0=0.5;
249  Double_t e_T0_eta=0.0;
250  Double_t e_n=2.0;
251  Double_t e_n_eta=0.0;
252  Double_t eta_p2=0.0;
253  Double_t etaMax=4.0;
254  if(isData) {
255  e_T0=0.6;
256  e_n=2.5;
257  e_T0_eta=0.05;
258  e_n_eta=-0.05;
259  eta_p2=1.5;
260  }
261  if(eta_p2>0.0) {
262  fEtaGen=TMath::Power(rnd->Uniform(),eta_p2)*etaMax;
263  if(rnd->Uniform()>=0.5) fEtaGen= -fEtaGen;
264  } else {
265  fEtaGen=rnd->Uniform(-etaMax,etaMax);
266  }
267  Double_t n=e_n + e_n_eta*fEtaGen;
268  Double_t T0=e_T0 + e_T0_eta*fEtaGen;
269  fPtGen=(TMath::Power(rnd->Uniform(),1./(1.-n))-1.)*T0;
270  /* static int print=100;
271  if(print) {
272  cout<<fEtaGen
273  <<" "<<fPtGen
274  <<"\n";
275  print--;
276  } */
277 }
278 
279 void ToyEvent::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
280  fPtGen=0.0;
281  fEtaGen=0.0;
282  fPtRec=rnd->Exp(2.5);
283  fEtaRec=rnd->Uniform(-3.,3.);
284 }
285 
286 void ToyEvent::GenerateReco(TRandom *rnd) {
287  if(fIsSignal) {
288  Double_t expEta=TMath::Exp(fEtaGen);
289  Double_t eGen=fPtGen*(expEta+1./expEta);
290  Double_t sigmaE=0.1*TMath::Sqrt(eGen)+(0.01+0.002*TMath::Abs(fEtaGen))
291  *eGen;
292  Double_t eRec;
293  do {
294  eRec=rnd->Gaus(eGen,sigmaE);
295  } while(eRec<=0.0);
296  Double_t sigmaEta=0.1+0.02*fEtaGen;
297  fEtaRec=rnd->Gaus(fEtaGen,sigmaEta);
298  fPtRec=eRec/(expEta+1./expEta);
299  do {
300  Double_t tauDiscr=0.08-0.04/(1.+fPtRec/10.0);
301  Double_t sigmaDiscr=0.01;
302  fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
303  } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
304  /* static int print=100;
305  if(print) {
306  cout<<fEtaGen<<" "<<fPtGen
307  <<" -> "<<fEtaRec<<" "<<fPtRec
308  <<"\n";
309  print--;
310  } */
311  } else {
312  do {
313  Double_t tauDiscr=0.15-0.05/(1.+fPtRec/5.0)+0.1*fEtaRec;
314  Double_t sigmaDiscr=0.02+0.01*fEtaRec;
315  fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
316  } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
317  }
318  fIsTriggered=(rnd->Uniform()<1./(TMath::Exp(-fPtRec+3.5)+1.));
319 }
Random number generator class based on M.
Definition: TRandom3.h:29
float Float_t
Definition: RtypesCore.h:53
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:235
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
T etaMax()
Function providing the maximum possible value of pseudorapidity for a non-zero rho, in the Scalar type with the largest dynamic range.
Definition: etaMax.h:50
const Bool_t kFALSE
Definition: Rtypes.h:92
STL namespace.
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:31
Double_t Exp(Double_t x)
Definition: TMath.h:495
double Double_t
Definition: RtypesCore.h:55
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
const Bool_t kTRUE
Definition: Rtypes.h:91
const Int_t n
Definition: legend1.C:16
virtual Double_t Exp(Double_t tau)
Returns an exponential deviate.
Definition: TRandom.cxx:212