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
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
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:22 PM.
/* 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 <cmath>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TDOMParser.h>
#include <TXMLDocument.h>
#include "TUnfoldBinningXML.h"
using std::cout;
switch on histogram errors
TH1::SetDefaultSumw2();
======================================================= 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;
TString dir = gSystem->UnixPathName(gSystem->GetDirName(__FILE__));
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();
TH2 *histMcsigGenRecF=TUnfoldBinning::CreateHistogramOfMigrations
(coarseBinning,fineBinning,"histMcsigGenRecF");
TH2 *histMcsigGenRecC=TUnfoldBinning::CreateHistogramOfMigrations
(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->ResetBranchAddresses();
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;
Draw all canvases
gROOT->GetListOfCanvases()->Draw()