Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
testUnfold7b.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 example is documented in conference proceedings:

arXiv:1611.01927 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)

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

The example comprises several macros

  • testUnfold7a.C create root files with TTree objects for signal, background and data
    • write files testUnfold7_signal.root testUnfold7_background.root testUnfold7_data.root
  • testUnfold7b.C loop over trees and fill histograms based on the TUnfoldBinning objects
    • read testUnfold7binning.xml testUnfold7_signal.root testUnfold7_background.root testUnfold7_data.root
    • write testUnfold7_histograms.root
  • testUnfold7c.C run the unfolding
    • read testUnfold7_histograms.root
    • write testUnfold7_result.root testUnfold7_result.ps
TUnfoldBinning "fine" has 33 bins [1,34] nTH1x=32
distribution: 33 bins
"pt" nbin=32 plus overflow
TUnfoldBinning "coarse" has 17 bins [1,18] nTH1x=16
distribution: 17 bins
"pt" nbin=16 plus overflow
loop over data events
loop over MC signal events
loop over MC background events
/* below is the content of the file testUnfold7binning.xml,
which is required as input to run this example.
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
<TUnfoldBinning>
<BinningNode name="fine" firstbin="1" factor="1">
<Axis name="pt" lowEdge="0.">
<Bin repeat="20" width="1."/>
<Bin repeat="12" width="2.5"/>
<Bin location="overflow" width="10"/>
</Axis>
</BinningNode>
<BinningNode name="coarse" firstbin="1" factor="1">
<Axis name="pt" lowEdge="0.">
<Bin repeat="10" width="2"/>
<Bin repeat="6" width="5"/>
<Bin location="overflow" width="10"/>
</Axis>
</BinningNode>
</TUnfoldBinning>
*/
#include <iostream>
#include <map>
#include <cmath>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TDOMParser.h>
#include <TXMLDocument.h>
using namespace std;
void testUnfold7b()
{
// switch on histogram errors
//=======================================================
// Step 1: open file to save histograms and binning schemes
TFile *outputFile=new TFile("testUnfold7_histograms.root","recreate");
//=======================================================
// Step 2: read binning from XML
// and save them to output file
TUnfoldBinning *fineBinningRoot,*coarseBinningRoot;
outputFile->cd();
// read binning schemes in XML format
TDOMParser parser;
Int_t error = parser.ParseFile(dir+"/testUnfold7binning.xml");
if(error) {
cout<<"error="<<error<<" from TDOMParser\n";
cout<<"==============================================================\n";
cout<<"Maybe the file testUnfold7binning.xml is missing?\n";
cout<<"The content of the file is included in the comments section\n";
cout<<"of this macro \"testUnfold7b.C\"\n";
cout<<"==============================================================\n";
}
TXMLDocument const *XMLdocument=parser.GetXMLDocument();
fineBinningRoot=
TUnfoldBinningXML::ImportXML(XMLdocument,"fine");
coarseBinningRoot=
TUnfoldBinningXML::ImportXML(XMLdocument,"coarse");
// write binning schemes to output file
fineBinningRoot->Write();
coarseBinningRoot->Write();
if(fineBinningRoot) {
fineBinningRoot->PrintStream(cout);
} else {
cout<<"could not read 'detector' binning\n";
}
if(coarseBinningRoot) {
coarseBinningRoot->PrintStream(cout);
} else {
cout<<"could not read 'generator' binning\n";
}
TUnfoldBinning const *fineBinning=fineBinningRoot;//->FindNode("ptfine");
TUnfoldBinning const *coarseBinning=coarseBinningRoot;//->FindNode("ptcoarse");
//=======================================================
// Step 3: book and fill data histograms
Float_t ptRec[3],ptGen[3],weight;
Int_t isTriggered,isSignal;
outputFile->cd();
TH1 *histDataRecF=fineBinning->CreateHistogram("histDataRecF");
TH1 *histDataRecC=coarseBinning->CreateHistogram("histDataRecC");
TH1 *histDataBgrF=fineBinning->CreateHistogram("histDataBgrF");
TH1 *histDataBgrC=coarseBinning->CreateHistogram("histDataBgrC");
TH1 *histDataGen=coarseBinning->CreateHistogram("histDataGen");
TFile *dataFile=new TFile("testUnfold7_data.root");
TTree *dataTree=(TTree *) dataFile->Get("data");
if(!dataTree) {
cout<<"could not read 'data' tree\n";
}
dataTree->ResetBranchAddresses();
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("ptgen",ptGen);
dataTree->SetBranchAddress("issignal",&isSignal);
dataTree->SetBranchStatus("*",1);
cout<<"loop over data events\n";
#define VAR_REC (ptRec[2])
#define VAR_GEN (ptGen[2])
for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
if(dataTree->GetEntry(ievent)<=0) break;
// fill histogram with reconstructed quantities
if(isTriggered) {
int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
histDataRecF->Fill(binF);
histDataRecC->Fill(binC);
if(!isSignal) {
histDataBgrF->Fill(binF);
histDataBgrC->Fill(binC);
}
}
// fill histogram with data truth parameters
if(isSignal) {
int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
histDataGen->Fill(binGen);
}
}
delete dataTree;
delete dataFile;
//=======================================================
// Step 4: book and fill histogram of migrations
// it receives events from both signal MC and background MC
outputFile->cd();
(coarseBinning,fineBinning,"histMcsigGenRecF");
(coarseBinning,coarseBinning,"histMcsigGenRecC");
TH1 *histMcsigRecF=fineBinning->CreateHistogram("histMcsigRecF");
TH1 *histMcsigRecC=coarseBinning->CreateHistogram("histMcsigRecC");
TH1 *histMcsigGen=coarseBinning->CreateHistogram("histMcsigGen");
TFile *signalFile=new TFile("testUnfold7_signal.root");
TTree *signalTree=(TTree *) signalFile->Get("signal");
if(!signalTree) {
cout<<"could not read 'signal' tree\n";
}
signalTree->ResetBranchAddresses();
signalTree->SetBranchAddress("ptrec",&ptRec);
//signalTree->SetBranchAddress("discr",&discr);
signalTree->SetBranchAddress("istriggered",&isTriggered);
signalTree->SetBranchAddress("ptgen",&ptGen);
signalTree->SetBranchAddress("weight",&weight);
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;
int binC=0,binF=0;
if(isTriggered) {
binF=fineBinning->GetGlobalBinNumber(VAR_REC);
binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
}
int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
histMcsigGenRecF->Fill(binGen,binF,weight);
histMcsigGenRecC->Fill(binGen,binC,weight);
histMcsigRecF->Fill(binF,weight);
histMcsigRecC->Fill(binC,weight);
histMcsigGen->Fill(binGen,weight);
}
delete signalTree;
delete signalFile;
outputFile->cd();
TH1 *histMcbgrRecF=fineBinning->CreateHistogram("histMcbgrRecF");
TH1 *histMcbgrRecC=coarseBinning->CreateHistogram("histMcbgrRecC");
TFile *bgrFile=new TFile("testUnfold7_background.root");
TTree *bgrTree=(TTree *) bgrFile->Get("background");
if(!bgrTree) {
cout<<"could not read 'background' tree\n";
}
bgrTree->SetBranchAddress("ptrec",&ptRec);
//bgrTree->SetBranchAddress("discr",&discr);
bgrTree->SetBranchAddress("istriggered",&isTriggered);
bgrTree->SetBranchAddress("weight",&weight);
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) {
int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
histMcbgrRecF->Fill(binF,weight);
histMcbgrRecC->Fill(binC,weight);
}
}
delete bgrTree;
delete bgrFile;
outputFile->Write();
delete outputFile;
}
int Int_t
Definition RtypesCore.h:45
float Float_t
Definition RtypesCore.h:57
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
virtual TXMLDocument * GetXMLDocument() const
Returns the TXMLDocument.
virtual Int_t ParseFile(const char *filename)
Parse the XML file where filename is the XML file name.
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:54
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:2352
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
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:6663
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:798
Basic string class.
Definition TString.h:136
virtual const char * UnixPathName(const char *unixpathname)
Convert from a local pathname to a Unix pathname.
Definition TSystem.cxx:1061
virtual TString GetDirName(const char *pathname)
Return the directory name in pathname.
Definition TSystem.cxx:1030
A TTree represents a columnar dataset.
Definition TTree.h:79
virtual Int_t SetBranchAddress(const char *bname, void *add, TBranch **ptr=0)
Change branch address, dealing with clone trees properly.
Definition TTree.cxx:8349
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:5618
virtual Long64_t GetEntriesFast() const
Definition TTree.h:462
virtual void ResetBranchAddresses()
Tell all of our branches to drop their current objects and allocate new ones.
Definition TTree.cxx:8047
virtual void SetBranchStatus(const char *bname, Bool_t status=1, UInt_t *found=0)
Set branch status to Process or DoNotProcess.
Definition TTree.cxx:8498
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 ...

Version 17.6, in parallel to changes in TUnfold

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 testUnfold7b.C.