Logo ROOT  
Reference Guide
testUnfold7b.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
60/* below is the content of the file testUnfold7binning.xml,
61 which is required as input to run this example.
62
63<?xml version="1.0" encoding="UTF-8" standalone="no"?>
64<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
65<TUnfoldBinning>
66<BinningNode name="fine" firstbin="1" factor="1">
67 <Axis name="pt" lowEdge="0.">
68 <Bin repeat="20" width="1."/>
69 <Bin repeat="12" width="2.5"/>
70 <Bin location="overflow" width="10"/>
71 </Axis>
72</BinningNode>
73<BinningNode name="coarse" firstbin="1" factor="1">
74 <Axis name="pt" lowEdge="0.">
75 <Bin repeat="10" width="2"/>
76 <Bin repeat="6" width="5"/>
77 <Bin location="overflow" width="10"/>
78 </Axis>
79</BinningNode>
80</TUnfoldBinning>
81
82 */
83
84#include <iostream>
85#include <map>
86#include <cmath>
87#include <TMath.h>
88#include <TFile.h>
89#include <TTree.h>
90#include <TH1.h>
91#include <TDOMParser.h>
92#include <TXMLDocument.h>
93#include "TUnfoldBinningXML.h"
94
95using namespace std;
96
97
98void testUnfold7b()
99{
100 // switch on histogram errors
102
103 //=======================================================
104 // Step 1: open file to save histograms and binning schemes
105
106 TFile *outputFile=new TFile("testUnfold7_histograms.root","recreate");
107
108 //=======================================================
109 // Step 2: read binning from XML
110 // and save them to output file
111
112 TUnfoldBinning *fineBinningRoot,*coarseBinningRoot;
113
114 outputFile->cd();
115
116 // read binning schemes in XML format
117
118 TDOMParser parser;
119 TString dir = gSystem->UnixPathName(gSystem->GetDirName(__FILE__));
120 Int_t error = parser.ParseFile(dir+"/testUnfold7binning.xml");
121 if(error) {
122 cout<<"error="<<error<<" from TDOMParser\n";
123 cout<<"==============================================================\n";
124 cout<<"Maybe the file testUnfold7binning.xml is missing?\n";
125 cout<<"The content of the file is included in the comments section\n";
126 cout<<"of this macro \"testUnfold7b.C\"\n";
127 cout<<"==============================================================\n";
128 }
129 TXMLDocument const *XMLdocument=parser.GetXMLDocument();
130 fineBinningRoot=
131 TUnfoldBinningXML::ImportXML(XMLdocument,"fine");
132 coarseBinningRoot=
133 TUnfoldBinningXML::ImportXML(XMLdocument,"coarse");
134
135 // write binning schemes to output file
136 fineBinningRoot->Write();
137 coarseBinningRoot->Write();
138
139 if(fineBinningRoot) {
140 fineBinningRoot->PrintStream(cout);
141 } else {
142 cout<<"could not read 'detector' binning\n";
143 }
144 if(coarseBinningRoot) {
145 coarseBinningRoot->PrintStream(cout);
146 } else {
147 cout<<"could not read 'generator' binning\n";
148 }
149
150 TUnfoldBinning const *fineBinning=fineBinningRoot;//->FindNode("ptfine");
151 TUnfoldBinning const *coarseBinning=coarseBinningRoot;//->FindNode("ptcoarse");
152
153
154 //=======================================================
155 // Step 3: book and fill data histograms
156
157 Float_t ptRec[3],ptGen[3],weight;
158 Int_t isTriggered,isSignal;
159
160 outputFile->cd();
161
162 TH1 *histDataRecF=fineBinning->CreateHistogram("histDataRecF");
163 TH1 *histDataRecC=coarseBinning->CreateHistogram("histDataRecC");
164 TH1 *histDataBgrF=fineBinning->CreateHistogram("histDataBgrF");
165 TH1 *histDataBgrC=coarseBinning->CreateHistogram("histDataBgrC");
166 TH1 *histDataGen=coarseBinning->CreateHistogram("histDataGen");
167
168 TFile *dataFile=new TFile("testUnfold7_data.root");
169 TTree *dataTree=(TTree *) dataFile->Get("data");
170
171 if(!dataTree) {
172 cout<<"could not read 'data' tree\n";
173 }
174
175 dataTree->ResetBranchAddresses();
176 dataTree->SetBranchAddress("ptrec",ptRec);
177 //dataTree->SetBranchAddress("discr",&discr);
178 // for real data, only the triggered events are available
179 dataTree->SetBranchAddress("istriggered",&isTriggered);
180 // data truth parameters
181 dataTree->SetBranchAddress("ptgen",ptGen);
182 dataTree->SetBranchAddress("issignal",&isSignal);
183 dataTree->SetBranchStatus("*",1);
184
185 cout<<"loop over data events\n";
186
187#define VAR_REC (ptRec[2])
188#define VAR_GEN (ptGen[2])
189
190 for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
191 if(dataTree->GetEntry(ievent)<=0) break;
192 // fill histogram with reconstructed quantities
193 if(isTriggered) {
194 int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
195 int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
196 histDataRecF->Fill(binF);
197 histDataRecC->Fill(binC);
198 if(!isSignal) {
199 histDataBgrF->Fill(binF);
200 histDataBgrC->Fill(binC);
201 }
202 }
203 // fill histogram with data truth parameters
204 if(isSignal) {
205 int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
206 histDataGen->Fill(binGen);
207 }
208 }
209
210 delete dataTree;
211 delete dataFile;
212
213 //=======================================================
214 // Step 4: book and fill histogram of migrations
215 // it receives events from both signal MC and background MC
216
217 outputFile->cd();
218
220 (coarseBinning,fineBinning,"histMcsigGenRecF");
222 (coarseBinning,coarseBinning,"histMcsigGenRecC");
223 TH1 *histMcsigRecF=fineBinning->CreateHistogram("histMcsigRecF");
224 TH1 *histMcsigRecC=coarseBinning->CreateHistogram("histMcsigRecC");
225 TH1 *histMcsigGen=coarseBinning->CreateHistogram("histMcsigGen");
226
227 TFile *signalFile=new TFile("testUnfold7_signal.root");
228 TTree *signalTree=(TTree *) signalFile->Get("signal");
229
230 if(!signalTree) {
231 cout<<"could not read 'signal' tree\n";
232 }
233
234 signalTree->ResetBranchAddresses();
235 signalTree->SetBranchAddress("ptrec",&ptRec);
236 //signalTree->SetBranchAddress("discr",&discr);
237 signalTree->SetBranchAddress("istriggered",&isTriggered);
238 signalTree->SetBranchAddress("ptgen",&ptGen);
239 signalTree->SetBranchAddress("weight",&weight);
240 signalTree->SetBranchStatus("*",1);
241
242 cout<<"loop over MC signal events\n";
243
244 for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
245 if(signalTree->GetEntry(ievent)<=0) break;
246
247 int binC=0,binF=0;
248 if(isTriggered) {
249 binF=fineBinning->GetGlobalBinNumber(VAR_REC);
250 binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
251 }
252 int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
253 histMcsigGenRecF->Fill(binGen,binF,weight);
254 histMcsigGenRecC->Fill(binGen,binC,weight);
255 histMcsigRecF->Fill(binF,weight);
256 histMcsigRecC->Fill(binC,weight);
257 histMcsigGen->Fill(binGen,weight);
258 }
259
260 delete signalTree;
261 delete signalFile;
262
263 outputFile->cd();
264
265 TH1 *histMcbgrRecF=fineBinning->CreateHistogram("histMcbgrRecF");
266 TH1 *histMcbgrRecC=coarseBinning->CreateHistogram("histMcbgrRecC");
267
268 TFile *bgrFile=new TFile("testUnfold7_background.root");
269 TTree *bgrTree=(TTree *) bgrFile->Get("background");
270
271 if(!bgrTree) {
272 cout<<"could not read 'background' tree\n";
273 }
274
275 bgrTree->ResetBranchAddresses();
276 bgrTree->SetBranchAddress("ptrec",&ptRec);
277 //bgrTree->SetBranchAddress("discr",&discr);
278 bgrTree->SetBranchAddress("istriggered",&isTriggered);
279 bgrTree->SetBranchAddress("weight",&weight);
280 bgrTree->SetBranchStatus("*",1);
281
282 cout<<"loop over MC background events\n";
283
284 for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
285 if(bgrTree->GetEntry(ievent)<=0) break;
286
287 // here, for background only reconstructed quantities are known
288 // and only the reconstructed events are relevant
289 if(isTriggered) {
290 int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
291 int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
292 histMcbgrRecF->Fill(binF,weight);
293 histMcbgrRecC->Fill(binC,weight);
294 }
295 }
296
297 delete bgrTree;
298 delete bgrFile;
299
300 outputFile->Write();
301 delete outputFile;
302
303}
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
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.
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
Basic string class.
Definition: TString.h:131
virtual const char * UnixPathName(const char *unixpathname)
Convert from a local pathname to a Unix pathname.
Definition: TSystem.cxx:1058
virtual TString GetDirName(const char *pathname)
Return the directory name in pathname.
Definition: TSystem.cxx:1027
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...
TXMLDocument contains a pointer to an xmlDoc structure, after the parser returns a tree built during ...
Definition: TXMLDocument.h:24