Logo ROOT  
Reference Guide
testUnfold5c.C File Reference

Detailed Description

View in nbviewer Open in SWAN Test program for the classes TUnfoldDensity and TUnfoldBinning.

A toy test of the TUnfold package

This is an example of unfolding a two-dimensional distribution also using an auxiliary measurement to constrain some background

The example comprises several macros

  • testUnfold5a.C create root files with TTree objects for signal, background and data
    • write files testUnfold5_signal.root testUnfold5_background.root testUnfold5_data.root
  • testUnfold5b.C create a root file with the TUnfoldBinning objects
    • write file testUnfold5_binning.root
  • testUnfold5c.C loop over trees and fill histograms based on the TUnfoldBinning objects
    • read testUnfold5_binning.root testUnfold5_signal.root testUnfold5_background.root testUnfold5_data.root
    • write testUnfold5_histograms.root
  • testUnfold5d.C run the unfolding
    • read testUnfold5_histograms.root
    • write testUnfold5_result.root testUnfold5_result.ps
TUnfoldBinning "detector" has 360 bins [1,361] nTH1x=360
TUnfoldBinning "detectordistribution" has 360 bins [1,361] nTH1x=360
distribution: 360 bins
"pt" nbin=8 plus overflow
"eta" nbin=10
"discriminator" nbin=4
TUnfoldBinning "generator" has 115 bins [1,116] nTH1x=115
TUnfoldBinning "signal" has 25 bins [1,26] nTH1x=25
distribution: 25 bins
"ptgen" nbin=3 plus underflow plus overflow
"etagen" nbin=3 plus underflow plus overflow
TUnfoldBinning "background" has 90 bins [26,116] nTH1x=90
distribution: 90 bins
"ptrec" nbin=8 plus overflow
"etarec" nbin=10
loop over data events
loop over MC signal events
loop over MC background events
// uncomment this to read the binning schemes from the root file
// by default the binning is read from the XML file
// #define READ_BINNING_CINT
#include <iostream>
#include <map>
#include <cmath>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#ifndef READ_BINNING_CINT
#include <TDOMParser.h>
#include <TXMLDocument.h>
#else
#include "TUnfoldBinning.h"
#endif
using namespace std;
void testUnfold5c()
{
// switch on histogram errors
//=======================================================
// Step 1: open file to save histograms and binning schemes
TFile *outputFile=new TFile("testUnfold5_histograms.root","recreate");
//=======================================================
// Step 2: read binning from XML
// and save them to output file
#ifdef READ_BINNING_CINT
TFile *binningSchemes=new TFile("testUnfold5_binning.root");
#endif
TUnfoldBinning *detectorBinning,*generatorBinning;
outputFile->cd();
// read binning schemes in XML format
#ifndef READ_BINNING_CINT
TDOMParser parser;
Int_t error=parser.ParseFile("testUnfold5binning.xml");
if(error) cout<<"error="<<error<<" from TDOMParser\n";
TXMLDocument const *XMLdocument=parser.GetXMLDocument();
detectorBinning=
TUnfoldBinningXML::ImportXML(XMLdocument,"detector");
generatorBinning=
TUnfoldBinningXML::ImportXML(XMLdocument,"generator");
#else
binningSchemes->GetObject("detector",detectorBinning);
binningSchemes->GetObject("generator",generatorBinning);
delete binningSchemes;
#endif
detectorBinning->Write();
generatorBinning->Write();
if(detectorBinning) {
detectorBinning->PrintStream(cout);
} else {
cout<<"could not read 'detector' binning\n";
}
if(generatorBinning) {
generatorBinning->PrintStream(cout);
} else {
cout<<"could not read 'generator' binning\n";
}
// pointers to various nodes in the binning scheme
const TUnfoldBinning *detectordistribution=
detectorBinning->FindNode("detectordistribution");
const TUnfoldBinning *signalBinning=
generatorBinning->FindNode("signal");
const TUnfoldBinning *bgrBinning=
generatorBinning->FindNode("background");
// write binning schemes to output file
//=======================================================
// Step 3: book and fill data histograms
Float_t etaRec,ptRec,discr,etaGen,ptGen;
Int_t istriggered,issignal;
outputFile->cd();
TH1 *histDataReco=detectorBinning->CreateHistogram("histDataReco");
TH1 *histDataTruth=generatorBinning->CreateHistogram("histDataTruth");
TFile *dataFile=new TFile("testUnfold5_data.root");
TTree *dataTree=(TTree *) dataFile->Get("data");
if(!dataTree) {
cout<<"could not read 'data' tree\n";
}
dataTree->ResetBranchAddresses();
dataTree->SetBranchAddress("etarec",&etaRec);
dataTree->SetBranchAddress("ptrec",&ptRec);
dataTree->SetBranchAddress("discr",&discr);
// for real data, only the triggered events are available
dataTree->SetBranchAddress("istriggered",&istriggered);
// data truth parameters
dataTree->SetBranchAddress("etagen",&etaGen);
dataTree->SetBranchAddress("ptgen",&ptGen);
dataTree->SetBranchAddress("issignal",&issignal);
dataTree->SetBranchStatus("*",1);
cout<<"loop over data events\n";
for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
if(dataTree->GetEntry(ievent)<=0) break;
// fill histogram with reconstructed quantities
if(istriggered) {
Int_t binNumber=
detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
histDataReco->Fill(binNumber);
}
// fill histogram with data truth parameters
if(issignal) {
// signal has true eta and pt
Int_t binNumber=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
histDataTruth->Fill(binNumber);
} else {
// background only has reconstructed pt and eta
Int_t binNumber=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
histDataTruth->Fill(binNumber);
}
}
delete dataTree;
delete dataFile;
//=======================================================
// Step 4: book and fill histogram of migrations
// it receives events from both signal MC and background MC
outputFile->cd();
(generatorBinning,detectorBinning,"histMCGenRec");
TFile *signalFile=new TFile("testUnfold5_signal.root");
TTree *signalTree=(TTree *) signalFile->Get("signal");
if(!signalTree) {
cout<<"could not read 'signal' tree\n";
}
signalTree->ResetBranchAddresses();
signalTree->SetBranchAddress("etarec",&etaRec);
signalTree->SetBranchAddress("ptrec",&ptRec);
signalTree->SetBranchAddress("discr",&discr);
signalTree->SetBranchAddress("istriggered",&istriggered);
signalTree->SetBranchAddress("etagen",&etaGen);
signalTree->SetBranchAddress("ptgen",&ptGen);
signalTree->SetBranchStatus("*",1);
cout<<"loop over MC signal events\n";
for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
if(signalTree->GetEntry(ievent)<=0) break;
// bin number on generator level for signal
Int_t genBin=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
// bin number on reconstructed level
// bin number 0 corresponds to non-reconstructed events
Int_t recBin=0;
if(istriggered) {
recBin=detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
}
histMCGenRec->Fill(genBin,recBin);
}
delete signalTree;
delete signalFile;
TFile *bgrFile=new TFile("testUnfold5_background.root");
TTree *bgrTree=(TTree *) bgrFile->Get("background");
if(!bgrTree) {
cout<<"could not read 'background' tree\n";
}
bgrTree->SetBranchAddress("etarec",&etaRec);
bgrTree->SetBranchAddress("ptrec",&ptRec);
bgrTree->SetBranchAddress("discr",&discr);
bgrTree->SetBranchAddress("istriggered",&istriggered);
bgrTree->SetBranchStatus("*",1);
cout<<"loop over MC background events\n";
for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
if(bgrTree->GetEntry(ievent)<=0) break;
// here, for background only reconstructed quantities are known
// and only the reconstructed events are relevant
if(istriggered) {
// bin number on generator level for background
Int_t genBin=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
// bin number on reconstructed level
Int_t recBin=detectordistribution->GetGlobalBinNumber
(ptRec,etaRec,discr);
histMCGenRec->Fill(genBin,recBin);
}
}
delete bgrTree;
delete bgrFile;
outputFile->Write();
delete outputFile;
}
int Int_t
Definition: RtypesCore.h:43
float Float_t
Definition: RtypesCore.h:55
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

Version 17.6, in parallel to changes in TUnfold

History:

  • Version 17.5, updated for reading binning from XML file
  • Version 17.4, updated for reading binning from XML file
  • Version 17.3, updated for reading binning from XML file
  • Version 17.2, updated for reading binning from XML file
  • Version 17.1, in parallel to changes in TUnfold
  • Version 17.0 example for multi-dimensional unfolding

This file is part of TUnfold.

TUnfold is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

TUnfold is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with TUnfold. If not, see http://www.gnu.org/licenses/.

Author
Stefan Schmitt DESY, 14.10.2008

Definition in file testUnfold5c.C.