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
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 Wednesday, April 17, 2024 at 11:24 AM.
%%cpp -d
#include <iostream>
#include <cmath>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TGraph.h>
#include <TFile.h>
#include <TH1.h>
#include "TUnfoldDensity.h"
using std::cout;
input_line_43:9:10: fatal error: 'TUnfoldDensity.h' file not found #include "TUnfoldDensity.h" ^~~~~~~~~~~~~~~~~~
switch on histogram errors
TH1::SetDefaultSumw2();
============================================== step 1 : open output file
TFile *outputFile=new TFile("testUnfold5_results.root","recreate");
============================================== step 2 : read binning schemes and input histograms
TFile *inputFile=new TFile("testUnfold5_histograms.root");
outputFile->cd();
TUnfoldBinning *detectorBinning,*generatorBinning;
inputFile->GetObject("detector",detectorBinning);
inputFile->GetObject("generator",generatorBinning);
if((!detectorBinning)||(!generatorBinning)) {
cout<<"problem to read binning schemes\n";
}
Error in <TFile::ReadBuffer>: error reading all requested bytes from file testUnfold5_histograms.root, got 250 of 300 Error in <TFile::Init>: testUnfold5_histograms.root failed to read the file type data. input_line_57:2:3: error: use of undeclared identifier 'TUnfoldBinning' (TUnfoldBinning * detectorBinning , * generatorBinning) ^ input_line_57:2:20: error: use of undeclared identifier 'detectorBinning' (TUnfoldBinning * detectorBinning , * generatorBinning) ^ input_line_57:2:40: error: use of undeclared identifier 'generatorBinning' (TUnfoldBinning * detectorBinning , * generatorBinning) ^ Error in <HandleInterpreterException>: Error evaluating expression (TUnfoldBinning * detectorBinning , * generatorBinning) Execution of your code was aborted.
save binning schemes to output file
detectorBinning->Write();
generatorBinning->Write();
input_line_59:2:3: error: use of undeclared identifier 'detectorBinning' (detectorBinning->Write()) ^ Error in <HandleInterpreterException>: Error evaluating expression (detectorBinning->Write()) Execution of your code was aborted.
read histograms
TH1 *histDataReco,*histDataTruth;
TH2 *histMCGenRec;
inputFile->GetObject("histDataReco",histDataReco);
inputFile->GetObject("histDataTruth",histDataTruth);
inputFile->GetObject("histMCGenRec",histMCGenRec);
#ifdef TEST_ZERO_UNCORR_ERROR
Unbalanced braces. This cell was not processed.
special test (bug in version 17.2 and below) set all errors in hisMCGenRec to zero -> program will crash
for(int i=0;i<=histMCGenRec->GetNbinsX()+1;i++) {
for(int j=0;j<=histMCGenRec->GetNbinsY()+1;j++) {
histMCGenRec->SetBinError(i,j,0.0);
}
}
#endif
histDataReco->Write();
histDataTruth->Write();
histMCGenRec->Write();
if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
cout<<"problem to read input histograms\n";
}
input_line_61:7:2: error: #endif without #if #endif ^ input_line_61:9:2: error: use of undeclared identifier 'histDataReco' histDataReco->Write(); ^ input_line_61:10:2: error: use of undeclared identifier 'histDataTruth' histDataTruth->Write(); ^ input_line_61:11:2: error: use of undeclared identifier 'histMCGenRec' histMCGenRec->Write(); ^ input_line_61:13:7: error: use of undeclared identifier 'histDataReco' if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) { ^ input_line_61:13:24: error: use of undeclared identifier 'histDataTruth' if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) { ^ input_line_61:13:42: error: use of undeclared identifier 'histMCGenRec' if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) { ^
======================== Step 3: unfolding
preserve the area
TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea;
input_line_62:2:3: error: 'TUnfold' is not a class, namespace, or enumeration TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea; ^ input_line_62:2:3: note: 'TUnfold' declared here input_line_62:2:40: error: use of undeclared identifier 'TUnfold' TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea; ^
basic choice of regularisation scheme: curvature (second derivative)
TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature;
input_line_63:2:3: error: 'TUnfold' is not a class, namespace, or enumeration TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature; ^ input_line_63:2:3: note: 'TUnfold' declared here input_line_63:2:31: error: use of undeclared identifier 'TUnfold' TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature; ^
density flags
TUnfoldDensity::EDensityMode densityFlags=
TUnfoldDensity::kDensityModeBinWidth;
input_line_64:2:3: error: 'TUnfoldDensity' is not a class, namespace, or enumeration TUnfoldDensity::EDensityMode densityFlags= ^ input_line_64:2:3: note: 'TUnfoldDensity' declared here input_line_64:3:4: error: use of undeclared identifier 'TUnfoldDensity' TUnfoldDensity::kDensityModeBinWidth; ^
detailed steering for regularisation
const char *REGULARISATION_DISTRIBUTION=nullptr;
const char *REGULARISATION_AXISSTEERING="*[B]";
set up matrix of migrations
TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz,
regMode,constraintMode,densityFlags,
generatorBinning,detectorBinning,
REGULARISATION_DISTRIBUTION,
REGULARISATION_AXISSTEERING);
input_line_66:2:3: error: unknown type name 'TUnfoldDensity' TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz, ^ input_line_66:2:38: error: use of undeclared identifier 'TUnfold' TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz, ^ input_line_66:2:25: error: use of undeclared identifier 'histMCGenRec' TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz, ^ input_line_66:3:24: error: use of undeclared identifier 'regMode' regMode,constraintMode,densityFlags, ^ input_line_66:3:32: error: use of undeclared identifier 'constraintMode' regMode,constraintMode,densityFlags, ^ input_line_66:3:47: error: use of undeclared identifier 'densityFlags' regMode,constraintMode,densityFlags, ^ input_line_66:4:24: error: use of undeclared identifier 'generatorBinning' generatorBinning,detectorBinning, ^ input_line_66:4:41: error: use of undeclared identifier 'detectorBinning' generatorBinning,detectorBinning, ^
define the input vector (the measured data distribution)
#ifdef TEST_INPUT_COVARIANCE
Unbalanced braces. This cell was not processed.
special test to use input covariance matrix
TH2D *inputEmatrix=
detectorBinning->CreateErrorMatrixHistogram("input_covar",true);
for(int i=1;i<=inputEmatrix->GetNbinsX();i++) {
Double_t e=histDataReco->GetBinError(i);
inputEmatrix->SetBinContent(i,i,e*e);
// test: non-zero covariance where variance is zero
// if(e<=0.) inputEmatrix->SetBinContent(i,i+1,1.0);
}
unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);
#else
unfold.SetInput(histDataReco /* ,0.0,1.0 */);
#endif
input_line_68:11:2: error: #else without #if #else ^ input_line_68:12:2: error: use of undeclared identifier 'unfold' unfold.SetInput(histDataReco /* ,0.0,1.0 */); ^ input_line_68:12:18: error: use of undeclared identifier 'histDataReco' unfold.SetInput(histDataReco /* ,0.0,1.0 */); ^ input_line_68:13:2: error: #endif without #if #endif ^
print matrix of regularisation conditions
#ifdef PRINT_MATRIX_L
TH2 *histL= unfold.GetL("L");
for(Int_t j=1;j<=histL->GetNbinsY();j++) {
cout<<"L["<<unfold.GetLBinning()->GetBinName(j)<<"]";
for(Int_t i=1;i<=histL->GetNbinsX();i++) {
Double_t c=histL->GetBinContent(i,j);
if(c!=0.0) cout<<" ["<<i<<"]="<<c;
}
cout<<"\n";
}
#endif
In module 'Cling_Runtime_Extra': /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:164:17: warning: inline function 'cling::Value::CastFwd<TH2D *>::cast' is not defined [-Wundefined-inline] static T* cast(const Value& V) { ^ /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:295:26: note: used here return CastFwd<T>::cast(*this); ^ /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:150:16: warning: inline function 'cling::Value::CastFwd<double>::cast' is not defined [-Wundefined-inline] static T cast(const Value& V) { ^ /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:295:26: note: used here return CastFwd<T>::cast(*this); ^
run the unfolding
here, tau is determined by scanning the global correlation coefficients
Int_t nScan=30;
TSpline *rhoLogTau=nullptr;
TGraph *lCurve=nullptr;
for determining tau, scan the correlation coefficients correlation coefficients may be probed for all distributions or only for selected distributions underflow/overflow bins may be included/excluded
const char *SCAN_DISTRIBUTION="signal";
const char *SCAN_AXISSTEERING=nullptr;
Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
TUnfoldDensity::kEScanTauRhoMax,
SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
&lCurve);
input_line_71:6:29: error: 'TUnfoldDensity' is not a class, namespace, or enumeration TUnfoldDensity::kEScanTauRhoMax, ^ input_line_71:6:29: note: 'TUnfoldDensity' declared here
create graphs with one point to visualize best choice of tau
Double_t t[1],rho[1],x[1],y[1];
rhoLogTau->GetKnot(iBest,t[0],rho[0]);
lCurve->GetPoint(iBest,x[0],y[0]);
TGraph *bestRhoLogTau=new TGraph(1,t,rho);
TGraph *bestLCurve=new TGraph(1,x,y);
Double_t *tAll=new Double_t[nScan],*rhoAll=new Double_t[nScan];
for(Int_t i=0;i<nScan;i++) {
rhoLogTau->GetKnot(i,tAll[i],rhoAll[i]);
}
TGraph *knots=new TGraph(nScan,tAll,rhoAll);
cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
<<" / "<<unfold.GetNdf()<<"\n";
input_line_73:2:44: error: use of undeclared identifier 'iBest' (((*(TSpline **)0x7fefc403a000))->GetKnot(iBest, (((Double_t(*))0x7fefc4037000))[0], (((Double_t(*))0x7fefc4037008))[0])) ^ Error in <HandleInterpreterException>: Error evaluating expression (((*(TSpline **)0x7fefc403a000))->GetKnot(iBest, (((Double_t(*))0x7fefc4037000))[0], (((Double_t(*))0x7fefc4037008))[0])) Execution of your code was aborted.
=========================== Step 4: retrieve and plot unfolding results
get unfolding output
TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",nullptr,nullptr,nullptr,kFALSE);
input_line_74:2:82: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Bool_t *' (aka 'const bool *') TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",nullptr,nullptr,nullptr,kFALSE); ^~~~~~
get Monte Carlo reconstructed data
TH1 *histMCReco=histMCGenRec->ProjectionY("histMCReco",0,-1,"e");
TH1 *histMCTruth=histMCGenRec->ProjectionX("histMCTruth",0,-1,"e");
Double_t scaleFactor=histDataTruth->GetSumOfWeights()/
histMCTruth->GetSumOfWeights();
histMCReco->Scale(scaleFactor);
histMCTruth->Scale(scaleFactor);
In module 'Cling_Runtime_Extra': /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:164:17: warning: inline function 'cling::Value::CastFwd<TH1 *>::cast' is not defined [-Wundefined-inline] static T* cast(const Value& V) { ^ /home/sftnight/build/workspace/root-makedoc-master/rootspi/rdoc/src/master.build/etc//cling/Interpreter/Value.h:295:26: note: used here return CastFwd<T>::cast(*this); ^ input_line_76:2:3: error: use of undeclared identifier 'histMCGenRec' (histMCGenRec->ProjectionY("histMCReco", 0, -1, "e")) ^ Error in <HandleInterpreterException>: Error evaluating expression (histMCGenRec->ProjectionY("histMCReco", 0, -1, "e")) Execution of your code was aborted.
get matrix of probabilities
TH2 *histProbability=unfold.GetProbabilityMatrix("histProbability");
input_line_78:2:3: error: use of undeclared identifier 'unfold' (unfold.GetProbabilityMatrix("histProbability")) ^ Error in <HandleInterpreterException>: Error evaluating expression (unfold.GetProbabilityMatrix("histProbability")) Execution of your code was aborted.
get global correlation coefficients
/* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",nullptr,nullptr,nullptr,kFALSE);
TH1 *histGlobalCorrScan=unfold.GetRhoItotal
("histGlobalCorrScan",nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);
/* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",nullptr,nullptr,nullptr,kFALSE);
TCanvas canvas;
canvas.Print("testUnfold5.ps[");
input_line_79:2:90: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Bool_t *' (aka 'const bool *') /* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",nullptr,nullptr,nullptr,kFALSE); ^~~~~~ input_line_79:4:71: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Bool_t *' (aka 'const bool *') ("histGlobalCorrScan",nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE); ^~~~~~ input_line_79:5:88: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Bool_t *' (aka 'const bool *') /* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",nullptr,nullptr,nullptr,kFALSE); ^~~~~~
========== page 1 ============ unfolding control plots input, matrix, output tau-scan, global correlations, correlation coefficients
canvas.Clear();
canvas.Divide(3,2);
input_line_81:2:3: error: use of undeclared identifier 'canvas' (canvas.Clear()) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.Clear()) Execution of your code was aborted.
(1) all bins, compare to original MC distribution
canvas.cd(1);
histDataReco->SetMinimum(0.0);
histDataReco->Draw("E");
histMCReco->SetLineColor(kBlue);
histMCReco->Draw("SAME HIST");
input_line_83:2:3: error: use of undeclared identifier 'canvas' (canvas.cd(1)) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.cd(1)) Execution of your code was aborted.
(2) matrix of probabilities
canvas.cd(2);
histProbability->Draw("BOX");
input_line_85:2:3: error: use of undeclared identifier 'canvas' (canvas.cd(2)) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.cd(2)) Execution of your code was aborted.
(3) unfolded data, data truth, MC truth
canvas.cd(3);
gPad->SetLogy();
histDataUnfold->Draw("E");
histDataTruth->SetLineColor(kBlue);
histDataTruth->Draw("SAME HIST");
histMCTruth->SetLineColor(kRed);
histMCTruth->Draw("SAME HIST");
input_line_86:5:30: error: cannot take the address of an rvalue of type 'EColor' histDataTruth->SetLineColor(kBlue); ^~~~~ Error while creating dynamic expression for: histDataTruth->SetLineColor(kBlue)
(4) scan of correlation vs tau
canvas.cd(4);
rhoLogTau->Draw();
knots->Draw("*");
bestRhoLogTau->SetMarkerColor(kRed);
bestRhoLogTau->Draw("*");
input_line_88:2:3: error: use of undeclared identifier 'canvas' (canvas.cd(4)) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.cd(4)) Execution of your code was aborted.
(5) global correlation coefficients for the distributions used during the scan
canvas.cd(5);
input_line_90:2:3: error: use of undeclared identifier 'canvas' (canvas.cd(5)) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.cd(5)) Execution of your code was aborted.
histCorrCoeff->Draw("BOX");
histGlobalCorrScan->Draw("HIST");
input_line_92:2:3: error: use of undeclared identifier 'histGlobalCorrScan' (histGlobalCorrScan->Draw("HIST")) ^ Error in <HandleInterpreterException>: Error evaluating expression (histGlobalCorrScan->Draw("HIST")) Execution of your code was aborted.
(6) L-curve
canvas.cd(6);
lCurve->Draw("AL");
bestLCurve->SetMarkerColor(kRed);
bestLCurve->Draw("*");
canvas.Print("testUnfold5.ps");
canvas.Print("testUnfold5.ps]");
input_line_94:2:3: error: use of undeclared identifier 'canvas' (canvas.cd(6)) ^ Error in <HandleInterpreterException>: Error evaluating expression (canvas.cd(6)) Execution of your code was aborted.
Draw all canvases
gROOT->GetListOfCanvases()->Draw()