Logo ROOT  
Reference Guide
testUnfold5d.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

chi**2=289.624+47.9755 / 213
#include <iostream>
#include <cmath>
#include <map>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TGraph.h>
#include <TFile.h>
#include <TH1.h>
#include "TUnfoldDensity.h"
using namespace std;
// #define PRINT_MATRIX_L
#define TEST_INPUT_COVARIANCE
void testUnfold5d()
{
// switch on histogram errors
//==============================================
// 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";
}
// save binning schemes to output file
detectorBinning->Write();
generatorBinning->Write();
// read histograms
TH1 *histDataReco,*histDataTruth;
TH2 *histMCGenRec;
inputFile->GetObject("histDataReco",histDataReco);
inputFile->GetObject("histDataTruth",histDataTruth);
inputFile->GetObject("histMCGenRec",histMCGenRec);
#ifdef TEST_ZERO_UNCORR_ERROR
// 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";
}
//========================
// Step 3: unfolding
// preserve the area
// basic choice of regularisation scheme:
// curvature (second derivative)
// density flags
// detailed steering for regularisation
const char *REGULARISATION_DISTRIBUTION=0;
const char *REGULARISATION_AXISSTEERING="*[B]";
// set up matrix of migrations
regMode,constraintMode,densityFlags,
generatorBinning,detectorBinning,
REGULARISATION_DISTRIBUTION,
REGULARISATION_AXISSTEERING);
// define the input vector (the measured data distribution)
#ifdef TEST_INPUT_COVARIANCE
// 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
// 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
// run the unfolding
//
// here, tau is determined by scanning the global correlation coefficients
Int_t nScan=30;
TSpline *rhoLogTau=0;
TGraph *lCurve=0;
// 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=0;
Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
&lCurve);
// 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";
//===========================
// Step 4: retrieve and plot unfolding results
// get unfolding output
TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",0,0,0,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);
// get matrix of probabilities
TH2 *histProbability=unfold.GetProbabilityMatrix("histProbability");
// get global correlation coefficients
/* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",0,0,0,kFALSE);
TH1 *histGlobalCorrScan=unfold.GetRhoItotal
("histGlobalCorrScan",0,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);
/* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",0,0,0,kFALSE);
TCanvas canvas;
canvas.Print("testUnfold5.ps[");
//========== page 1 ============
// unfolding control plots
// input, matrix, output
// tau-scan, global correlations, correlation coefficients
canvas.Clear();
canvas.Divide(3,2);
// (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");
// (2) matrix of probabilities
canvas.cd(2);
histProbability->Draw("BOX");
// (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");
// (4) scan of correlation vs tau
canvas.cd(4);
rhoLogTau->Draw();
knots->Draw("*");
bestRhoLogTau->SetMarkerColor(kRed);
bestRhoLogTau->Draw("*");
// (5) global correlation coefficients for the distributions
// used during the scan
canvas.cd(5);
//histCorrCoeff->Draw("BOX");
histGlobalCorrScan->Draw("HIST");
// (6) L-curve
canvas.cd(6);
lCurve->Draw("AL");
bestLCurve->SetMarkerColor(kRed);
bestLCurve->Draw("*");
canvas.Print("testUnfold5.ps");
canvas.Print("testUnfold5.ps]");
}
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
#define gPad
Definition: TVirtualPad.h:287
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
The Canvas class.
Definition: TCanvas.h:27
void Clear(Option_t *option="")
Remove all primitives from the canvas.
Definition: TCanvas.cxx:720
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:701
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
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
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:760
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
Definition: TGraph.cxx:1593
The TH1 histogram class.
Definition: TH1.h:56
virtual Int_t GetNbinsY() const
Definition: TH1.h:293
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8662
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:395
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
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6246
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7426
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along Y.
Definition: TH2.cxx:2340
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along X.
Definition: TH2.cxx:2300
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:88
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2480
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
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1165
virtual void Print(const char *filename="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:4667
Base class for spline implementation containing the Draw/Paint methods.
Definition: TSpline.h:22
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TSpline.cxx:97
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
Create a TH2D histogram capable to hold a covariance matrix.
An algorithm to unfold distributions from detector to truth level.
@ kEScanTauRhoMax
maximum global correlation coefficient (from TUnfold::GetRhoI())
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeBinWidth
scale factors from multidimensional bin width
EConstraint
type of extra constraint
Definition: TUnfold.h:109
@ kEConstraintArea
enforce preservation of the area
Definition: TUnfold.h:115
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
Definition: TUnfold.h:131
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17

Version 17.6, in parallel to changes in TUnfold

History:

  • Version 17.5, in parallel to changes in TUnfold
  • Version 17.4, in parallel to changes in TUnfold
  • Version 17.3, in parallel to changes in TUnfold
  • Version 17.2, in parallel to changes in TUnfold
  • 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 testUnfold5d.C.