Logo ROOT  
Reference Guide
testUnfold5c.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, updated for reading binning from XML file
42/// - Version 17.4, updated for reading binning from XML file
43/// - Version 17.3, updated for reading binning from XML file
44/// - Version 17.2, updated for reading binning from XML file
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// uncomment this to read the binning schemes from the root file
66// by default the binning is read from the XML file
67// #define READ_BINNING_CINT
68
69
70#include <iostream>
71#include <map>
72#include <cmath>
73#include <TMath.h>
74#include <TFile.h>
75#include <TTree.h>
76#include <TH1.h>
77#ifndef READ_BINNING_CINT
78#include <TDOMParser.h>
79#include <TXMLDocument.h>
80#include "TUnfoldBinningXML.h"
81#else
82#include "TUnfoldBinning.h"
83#endif
84
85using namespace std;
86
87
88void testUnfold5c()
89{
90 // switch on histogram errors
92
93 //=======================================================
94 // Step 1: open file to save histograms and binning schemes
95
96 TFile *outputFile=new TFile("testUnfold5_histograms.root","recreate");
97
98 //=======================================================
99 // Step 2: read binning from XML
100 // and save them to output file
101
102#ifdef READ_BINNING_CINT
103 TFile *binningSchemes=new TFile("testUnfold5_binning.root");
104#endif
105
106 TUnfoldBinning *detectorBinning,*generatorBinning;
107
108 outputFile->cd();
109
110 // read binning schemes in XML format
111#ifndef READ_BINNING_CINT
112 TDOMParser parser;
113 Int_t error=parser.ParseFile("testUnfold5binning.xml");
114 if(error) cout<<"error="<<error<<" from TDOMParser\n";
115 TXMLDocument const *XMLdocument=parser.GetXMLDocument();
116 detectorBinning=
117 TUnfoldBinningXML::ImportXML(XMLdocument,"detector");
118 generatorBinning=
119 TUnfoldBinningXML::ImportXML(XMLdocument,"generator");
120#else
121 binningSchemes->GetObject("detector",detectorBinning);
122 binningSchemes->GetObject("generator",generatorBinning);
123
124 delete binningSchemes;
125#endif
126 detectorBinning->Write();
127 generatorBinning->Write();
128
129 if(detectorBinning) {
130 detectorBinning->PrintStream(cout);
131 } else {
132 cout<<"could not read 'detector' binning\n";
133 }
134 if(generatorBinning) {
135 generatorBinning->PrintStream(cout);
136 } else {
137 cout<<"could not read 'generator' binning\n";
138 }
139
140 // pointers to various nodes in the binning scheme
141 const TUnfoldBinning *detectordistribution=
142 detectorBinning->FindNode("detectordistribution");
143
144 const TUnfoldBinning *signalBinning=
145 generatorBinning->FindNode("signal");
146
147 const TUnfoldBinning *bgrBinning=
148 generatorBinning->FindNode("background");
149
150 // write binning schemes to output file
151
152 //=======================================================
153 // Step 3: book and fill data histograms
154
155 Float_t etaRec,ptRec,discr,etaGen,ptGen;
156 Int_t istriggered,issignal;
157
158 outputFile->cd();
159
160 TH1 *histDataReco=detectorBinning->CreateHistogram("histDataReco");
161 TH1 *histDataTruth=generatorBinning->CreateHistogram("histDataTruth");
162
163 TFile *dataFile=new TFile("testUnfold5_data.root");
164 TTree *dataTree=(TTree *) dataFile->Get("data");
165
166 if(!dataTree) {
167 cout<<"could not read 'data' tree\n";
168 }
169
170 dataTree->ResetBranchAddresses();
171 dataTree->SetBranchAddress("etarec",&etaRec);
172 dataTree->SetBranchAddress("ptrec",&ptRec);
173 dataTree->SetBranchAddress("discr",&discr);
174 // for real data, only the triggered events are available
175 dataTree->SetBranchAddress("istriggered",&istriggered);
176 // data truth parameters
177 dataTree->SetBranchAddress("etagen",&etaGen);
178 dataTree->SetBranchAddress("ptgen",&ptGen);
179 dataTree->SetBranchAddress("issignal",&issignal);
180 dataTree->SetBranchStatus("*",1);
181
182
183 cout<<"loop over data events\n";
184
185 for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
186 if(dataTree->GetEntry(ievent)<=0) break;
187 // fill histogram with reconstructed quantities
188 if(istriggered) {
189 Int_t binNumber=
190 detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
191 histDataReco->Fill(binNumber);
192 }
193 // fill histogram with data truth parameters
194 if(issignal) {
195 // signal has true eta and pt
196 Int_t binNumber=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
197 histDataTruth->Fill(binNumber);
198 } else {
199 // background only has reconstructed pt and eta
200 Int_t binNumber=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
201 histDataTruth->Fill(binNumber);
202 }
203 }
204
205 delete dataTree;
206 delete dataFile;
207
208 //=======================================================
209 // Step 4: book and fill histogram of migrations
210 // it receives events from both signal MC and background MC
211
212 outputFile->cd();
213
215 (generatorBinning,detectorBinning,"histMCGenRec");
216
217 TFile *signalFile=new TFile("testUnfold5_signal.root");
218 TTree *signalTree=(TTree *) signalFile->Get("signal");
219
220 if(!signalTree) {
221 cout<<"could not read 'signal' tree\n";
222 }
223
224 signalTree->ResetBranchAddresses();
225 signalTree->SetBranchAddress("etarec",&etaRec);
226 signalTree->SetBranchAddress("ptrec",&ptRec);
227 signalTree->SetBranchAddress("discr",&discr);
228 signalTree->SetBranchAddress("istriggered",&istriggered);
229 signalTree->SetBranchAddress("etagen",&etaGen);
230 signalTree->SetBranchAddress("ptgen",&ptGen);
231 signalTree->SetBranchStatus("*",1);
232
233 cout<<"loop over MC signal events\n";
234
235 for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
236 if(signalTree->GetEntry(ievent)<=0) break;
237
238 // bin number on generator level for signal
239 Int_t genBin=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
240
241 // bin number on reconstructed level
242 // bin number 0 corresponds to non-reconstructed events
243 Int_t recBin=0;
244 if(istriggered) {
245 recBin=detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
246 }
247 histMCGenRec->Fill(genBin,recBin);
248 }
249
250 delete signalTree;
251 delete signalFile;
252
253 TFile *bgrFile=new TFile("testUnfold5_background.root");
254 TTree *bgrTree=(TTree *) bgrFile->Get("background");
255
256 if(!bgrTree) {
257 cout<<"could not read 'background' tree\n";
258 }
259
260 bgrTree->ResetBranchAddresses();
261 bgrTree->SetBranchAddress("etarec",&etaRec);
262 bgrTree->SetBranchAddress("ptrec",&ptRec);
263 bgrTree->SetBranchAddress("discr",&discr);
264 bgrTree->SetBranchAddress("istriggered",&istriggered);
265 bgrTree->SetBranchStatus("*",1);
266
267 cout<<"loop over MC background events\n";
268
269 for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
270 if(bgrTree->GetEntry(ievent)<=0) break;
271
272 // here, for background only reconstructed quantities are known
273 // and only the reconstructed events are relevant
274 if(istriggered) {
275 // bin number on generator level for background
276 Int_t genBin=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
277 // bin number on reconstructed level
278 Int_t recBin=detectordistribution->GetGlobalBinNumber
279 (ptRec,etaRec,discr);
280 histMCGenRec->Fill(genBin,recBin);
281 }
282 }
283
284 delete bgrTree;
285 delete bgrFile;
286
287 outputFile->Write();
288 delete outputFile;
289
290}
virtual TXMLDocument * GetXMLDocument() const
Returns the TXMLDocument.
Definition: TDOMParser.cxx:144
virtual Int_t ParseFile(const char *filename)
Parse the XML file where filename is the XML file name.
Definition: TDOMParser.cxx:70
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
void GetObject(const char *namecycle, T *&ptr)
Definition: TDirectory.h:155
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
Int_t Write(const char *name=nullptr, Int_t opt=0, Int_t bufsiz=0) override
Write memory objects to this file.
Definition: TFile.cxx:2297
The TH1 histogram class.
Definition: TH1.h:56
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
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:6330
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:796
A TTree represents a columnar dataset.
Definition: TTree.h:78
virtual Int_t SetBranchAddress(const char *bname, void *add, TBranch **ptr=0)
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:8237
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5542
virtual Long64_t GetEntriesFast() const
Definition: TTree.h:459
virtual void ResetBranchAddresses()
Tell all of our branches to drop their current objects and allocate new ones.
Definition: TTree.cxx:7969
virtual void SetBranchStatus(const char *bname, Bool_t status=1, UInt_t *found=0)
Set branch status to Process or DoNotProcess.
Definition: TTree.cxx:8386
static TUnfoldBinningXML * ImportXML(const TXMLDocument *document, const char *name)
Import a binning scheme from an XML file.
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
void PrintStream(std::ostream &out, Int_t indent=0, int debug=0) const
Print some information about this binning tree.
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
Create a THxx histogram capable to hold the bins of this binning node and its children.
Int_t GetGlobalBinNumber(Double_t x) const
Locate a bin in a one-dimensional distribution.
static TH2D * CreateHistogramOfMigrations(TUnfoldBinning const *xAxis, TUnfoldBinning const *yAxis, char const *histogramName, Bool_t originalXAxisBinning=kFALSE, Bool_t originalYAxisBinning=kFALSE, char const *histogramTitle=0)
Create a TH2D histogram capable to hold the bins of the two input binning schemes on the x and y axes...
TUnfoldBinning const * FindNode(char const *name) const
Traverse the tree and return the first node which matches the given name.
TXMLDocument contains a pointer to an xmlDoc structure, after the parser returns a tree built during ...
Definition: TXMLDocument.h:24