Logo ROOT  
Reference Guide
No Matches
testUnfold3.C File Reference

Detailed Description

View in nbviewer Open in SWAN Simple Test program for the class TUnfoldDensity.

1-dimensional unfolding with background subtraction

the goal is to unfold the underlying "true" distribution of a variable Pt

the reconstructed Pt is measured in 24 bins from 4 to 28 the generator-level Pt is unfolded into 10 bins from 6 to 26

  • plus underflow bin from 0 to 6
  • plus overflow bin above 26 there are two background sources bgr1 and bgr2 the signal has a finite trigger efficiency at a threshold of 8 GeV

one type of systematic error is studied, where the signal parameters are changed

Finally, the unfolding is compared to a "bin-by-bin" correction method

TUnfold version is V17.6
chi**2=25.3912+4.28619 / 11
bin truth result (stat) (bgr) (sys)
1 4002 4261.3 +/-137.0 +/- 8.9 +/-473.2 (unfolding)
4190.1 +/-128.7 +/-104.7 +/-266.3 (bin-by-bin)
2 1816 1775.1 +/- 82.0 +/- 7.3 +/- 25.2 (unfolding)
1714.1 +/- 61.0 +/- 17.0 +/- 70.8 (bin-by-bin)
3 1011 1020.9 +/- 60.1 +/- 6.7 +/- 24.8 (unfolding)
959.5 +/- 40.9 +/- 9.2 +/- 43.3 (bin-by-bin)
4 593 526.9 +/- 45.2 +/- 6.1 +/- 15.0 (unfolding)
506.2 +/- 27.9 +/- 6.6 +/- 45.2 (bin-by-bin)
5 379 371.6 +/- 38.6 +/- 6.5 +/- 19.8 (unfolding)
351.8 +/- 22.6 +/- 6.0 +/- 29.2 (bin-by-bin)
6 256 341.5 +/- 34.4 +/- 6.0 +/- 9.2 (unfolding)
282.5 +/- 19.2 +/- 5.2 +/- 46.6 (bin-by-bin)
7 193 182.0 +/- 30.1 +/- 5.7 +/- 7.5 (unfolding)
177.7 +/- 15.8 +/- 4.9 +/- 22.5 (bin-by-bin)
8 137 147.1 +/- 28.4 +/- 5.6 +/- 7.0 (unfolding)
133.6 +/- 13.9 +/- 4.5 +/- 20.5 (bin-by-bin)
9 123 109.6 +/- 26.3 +/- 5.4 +/- 10.4 (unfolding)
98.6 +/- 12.1 +/- 4.1 +/- 20.7 (bin-by-bin)
10 78 100.7 +/- 24.7 +/- 6.1 +/- 8.4 (unfolding)
89.7 +/- 13.2 +/- 5.0 +/- 17.1 (bin-by-bin)
#include <TMath.h>
#include <TCanvas.h>
#include <TRandom3.h>
#include <TFitter.h>
#include <TF1.h>
#include <TStyle.h>
#include <TVector.h>
#include <TGraph.h>
#include "TUnfoldDensity.h"
using namespace std;
TRandom *rnd=0;
Double_t GenerateEvent(const Double_t *parm,
const Double_t *triggerParm,
Double_t *intLumi,
Bool_t *triggerFlag,
Double_t *ptGen,Int_t *iType)
// generate an event
// input:
// parameters for the event generator
// return value:
// reconstructed Pt
// output to pointers:
// integrated luminosity
// several variables only accessible on generator level
// the parm array defines the physical parameters
// parm[0]: background source 1 fraction
// parm[1]: background source 2 fraction
// parm[2]: lower limit of generated Pt distribution
// parm[3]: upper limit of generated Pt distribution
// parm[4]: exponent for Pt distribution signal
// parm[5]: exponent for Pt distribution background 1
// parm[6]: exponent for Pt distribution background 2
// parm[7]: resolution parameter a goes with sqrt(Pt)
// parm[8]: resolution parameter b goes with Pt
// triggerParm[0]: trigger threshold turn-on position
// triggerParm[1]: trigger threshold turn-on width
// triggerParm[2]: trigger efficiency for high pt
// intLumi is advanced bu 1 for each *generated* event
// for data, several events may be generated, until one passes the trigger
// some generator-level quantities are also returned:
// triggerFlag: whether the event passed the trigger threshold
// ptGen: the generated pt
// iType: which type of process was simulated
// the "triggerFlag" also has another meaning:
// if(triggerFlag==0) only triggered events are returned
// Usage to generate data events
// ptObs=GenerateEvent(parm,triggerParm,0,0,0)
// Usage to generate MC events
// ptGen=GenerateEvent(parm,triggerParm,&triggerFlag,&ptGen,&iType);
Double_t ptObs;
Bool_t isTriggered=kFALSE;
do {
Int_t itype;
Double_t ptgen;
Double_t f=rnd->Rndm();
// decide whether this is background or signal
if(f<parm[0]) itype=1;
else if(f<parm[0]+parm[1]) itype=2;
// generate Pt according to distribution pow(ptgen,a)
// get exponent
Double_t a=parm[4+itype];
Double_t a1=a+1.0;
Double_t t=rnd->Rndm();
if(a1 == 0.0) {
Double_t x0=TMath::Log(parm[2]);
} else {
Double_t x0=pow(parm[2],a1);
if(iType) *iType=itype;
if(ptGen) *ptGen=ptgen;
// smearing in Pt with large asymmetric tail
// decide whether event was triggered
Double_t triggerProb =
isTriggered= rnd->Rndm()<triggerProb;
(*intLumi) ++;
} while((!triggerFlag) && (!isTriggered));
// return trigger decision
if(triggerFlag) *triggerFlag=isTriggered;
return ptObs;
//int main(int argc, char *argv[])
void testUnfold3()
// switch on histogram errors
// show fit result
// random generator
rnd=new TRandom3();
// data and MC luminosities
Double_t const lumiData= 30000;
Double_t const lumiMC =1000000;
// Step 1: define binning, book histograms
// reconstructed pt (fine binning)
Int_t const nDet=24;
Double_t const xminDet=4.0;
Double_t const xmaxDet=28.0;
// generated pt (coarse binning)
Int_t const nGen=10;
Double_t const xminGen= 6.0;
Double_t const xmaxGen=26.0;
// book histograms
// (1) unfolding input: binning scheme is fine for detector,coarse for gen
// histUnfoldInput : reconstructed data, binning for unfolding
// histUnfoldMatrix : generated vs reconstructed distribution
// histUnfoldBgr1 : background source1, as predicted by MC
// histUnfoldBgr2 : background source2, as predicted by MC
// for systematic studies
// histUnfoldMatrixSys : generated vs reconstructed with different signal
// (2) histograms required for bin-by-bin method
// histDetDATAbbb : reconstructed data for bin-by-bin
// histDetMCbbb : reconstructed MC for bin-by-bin
// histDetMCBgr1bbb : reconstructed MC bgr1 for bin-by-bin
// histDetMCBgr2bbb : reconstructed MC bgr2 for bin-by-bin
// histDetMCBgrPtbbb : reconstructed MC bgr from low/high pt for bin-by-bin
// histGenMCbbb : generated MC truth
// for systematic studies
// histDetMCSysbbb : reconstructed with changed trigger
// histGenMCSysbbb : generated MC truth
// (3) monitoring and control
// histGenData : data truth for bias tests
// histDetMC : MC prediction
// (1) create histograms required for unfolding
TH1D *histUnfoldInput=
new TH1D("unfolding input rec",";ptrec",nDet,xminDet,xmaxDet);
TH2D *histUnfoldMatrix=
new TH2D("unfolding matrix",";ptgen;ptrec",
TH1D *histUnfoldBgr1=
new TH1D("unfolding bgr1 rec",";ptrec",nDet,xminDet,xmaxDet);
TH1D *histUnfoldBgr2=
new TH1D("unfolding bgr2 rec",";ptrec",nDet,xminDet,xmaxDet);
TH2D *histUnfoldMatrixSys=
new TH2D("unfolding matrix sys",";ptgen;ptrec",
// (2) histograms required for bin-by-bin
TH1D *histBbbInput=
new TH1D("bbb input rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalRec=
new TH1D("bbb signal rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbBgr1=
new TH1D("bbb bgr1 rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbBgr2=
new TH1D("bbb bgr2 rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbBgrPt=
new TH1D("bbb bgrPt rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalGen=
new TH1D("bbb signal gen",";ptgen",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalRecSys=
new TH1D("bbb signal sys rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbBgrPtSys=
new TH1D("bbb bgrPt sys rec",";ptrec",nGen,xminGen,xmaxGen);
TH1D *histBbbSignalGenSys=
new TH1D("bbb signal sys gen",";ptgen",nGen,xminGen,xmaxGen);
// (3) control histograms
TH1D *histDataTruth=
new TH1D("DATA truth gen",";ptgen",nGen,xminGen,xmaxGen);
TH1D *histDetMC=
new TH1D("MC prediction rec",";ptrec",nDet,xminDet,xmaxDet);
// ==============================================================
// Step 2: generate data distributions
// data parameters: in real life these are unknown
static Double_t parm_DATA[]={
0.05, // fraction of background 1 (on generator level)
0.05, // fraction of background 2 (on generator level)
3.5, // lower Pt cut (generator level)
100.,// upper Pt cut (generator level)
-3.0,// signal exponent
0.1, // background 1 exponent
-0.5, // background 2 exponent
0.2, // energy resolution a term
0.01, // energy resolution b term
static Double_t triggerParm_DATA[]={8.0, // threshold is 8 GeV
4.0, // width is 4 GeV
0.95 // high Pt efficiency os 95%
Double_t intLumi=0.0;
while(intLumi<lumiData) {
Int_t iTypeGen;
Bool_t isTriggered;
Double_t ptGen;
Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
if(isTriggered) {
// (1) histogram for unfolding
// (2) histogram for bin-by-bin
// (3) monitoring
if(iTypeGen==0) histDataTruth->Fill(ptGen);
// ==============================================================
// Step 3: generate default MC distributions
// MC parameters
// default settings
static Double_t parm_MC[]={
0.05, // fraction of background 1 (on generator level)
0.05, // fraction of background 2 (on generator level)
3.5, // lower Pt cut (generator level)
100.,// upper Pt cut (generator level)
-4.0,// signal exponent !!! steeper than in data
// to illustrate bin-by-bin bias
0.1, // background 1 exponent
-0.5, // background 2 exponent
0.2, // energy resolution a term
0.01, // energy resolution b term
static Double_t triggerParm_MC[]={8.0, // threshold is 8 GeV
4.0, // width is 4 GeV
0.95 // high Pt efficiency is 95%
// weighting factor to accomodate for the different luminosity in data and MC
Double_t lumiWeight = lumiData/lumiMC;
while(intLumi<lumiMC) {
Int_t iTypeGen;
Bool_t isTriggered;
Double_t ptGen;
Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
if(!isTriggered) ptObs=0.0;
// (1) distribution required for unfolding
if(iTypeGen==0) {
} else if(iTypeGen==1) {
} else if(iTypeGen==2) {
// (2) distributions for bbb
if(iTypeGen==0) {
if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
} else {
} else if(iTypeGen==1) {
} else if(iTypeGen==2) {
// (3) control distribution
histDetMC ->Fill(ptObs,lumiWeight);
// ==============================================================
// Step 4: generate MC distributions for systematic study
// example: study dependence on initial signal shape
// -> BGR distributions do not change
static Double_t parm_MC_SYS[]={
0.05, // fraction of background: unchanged
0.05, // fraction of background: unchanged
3.5, // lower Pt cut (generator level)
100.,// upper Pt cut (generator level)
-2.0, // signal exponent changed
0.1, // background 1 exponent
-0.5, // background 2 exponent
0.2, // energy resolution a term
0.01, // energy resolution b term
while(intLumi<lumiMC) {
Int_t iTypeGen;
Bool_t isTriggered;
Double_t ptGen;
Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
if(!isTriggered) ptObs=0.0;
// (1) for unfolding
if(iTypeGen==0) {
// (2) for bin-by-bin
if(iTypeGen==0) {
if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
} else {
// this method is new in version 16 of TUnfold
cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
// Step 5: unfolding
// here one has to decide about the regularisation scheme
// regularize curvature
// preserve the area
TUnfold::EConstraint constraintMode=
// bin content is divided by the bin width
// set up matrix of migrations
// define the input vector (the measured data distribution)
// subtract background, normalized to data luminosity
// with 10% scale error each
Double_t scale_bgr=1.0;
Double_t dscale_bgr=0.1;
// add systematic error
// run the unfolding
Int_t nScan=30;
TSpline *logTauX,*logTauY;
TGraph *lCurve;
// this method scans the parameter tau and finds the kink in the L curve
// finally, the unfolding is done for the best choice of tau
Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
<<" / "<<unfold.GetNdf()<<"\n";
// save graphs with one point to visualize best choice of tau
Double_t t[1],x[1],y[1];
TGraph *bestLcurve=new TGraph(1,x,y);
TGraph *bestLogTauLogChi2=new TGraph(1,t,x);
// Step 6: retrieve unfolding results
// get unfolding output
// includes the statistical and background errors
// but not the other systematic uncertainties
TH1 *histUnfoldOutput=unfold.GetOutput("PT(unfold,stat+bgrerr)");
// retrieve error matrix of statistical errors
TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");
// retrieve full error matrix
// This includes all systematic errors
TH2 *histEmatTotal=unfold.GetEmatrixTotal("unfolding total error matrix");
// create two copies of the unfolded data, one with statistical errors
// the other with total errors
TH1 *histUnfoldStat=new TH1D("PT(unfold,staterr)",";Pt(gen)",
TH1 *histUnfoldTotal=new TH1D("PT(unfold,totalerr)",";Pt(gen)",
for(Int_t i=0;i<nGen+2;i++) {
Double_t c=histUnfoldOutput->GetBinContent(i);
// histogram with unfolded data and stat errors
// histogram with unfolded data and total errors
// create histogram with correlation matrix
TH2D *histCorr=new TH2D("Corr(total)",";Pt(gen);Pt(gen)",
for(Int_t i=0;i<nGen+2;i++) {
Double_t ei,ej;
if(ei<=0.0) continue;
for(Int_t j=0;j<nGen+2;j++) {
if(ej<=0.0) continue;
// retrieve bgr source 1
TH1 *histDetNormBgr1=unfold.GetBackground("bgr1 normalized",
TH1 *histDetNormBgrTotal=unfold.GetBackground("bgr total normalized");
// Step 7: plots
// data, MC prediction, background
histDetMC->Draw("SAME HIST");
histDetNormBgr1->Draw("SAME HIST");
histDetNormBgrTotal->Draw("SAME HIST");
// unfolded data, data truth, MC truth
// outer error: total error
// middle error: stat+bgr
histUnfoldOutput->Draw("SAME E1");
// inner error: stat only
histUnfoldStat->Draw("SAME E1");
histDataTruth->Draw("SAME HIST");
histBbbSignalGen->Draw("SAME HIST");
// unfolding matrix
// show tau as a function of chi**2
// show the L curve
// show correlation matrix
// step 8: compare results to the so-called bin-by-bin "correction"
std::cout<<"bin truth result (stat) (bgr) (sys)\n";
for(Int_t i=1;i<=nGen;i++) {
// data contribution in this bin
Double_t data=histBbbInput->GetBinContent(i);
Double_t errData=histBbbInput->GetBinError(i);
// subtract background contributions
Double_t data_bgr=data
- scale_bgr*(histBbbBgr1->GetBinContent(i)
+ histBbbBgr2->GetBinContent(i)
+ histBbbBgrPt->GetBinContent(i));
Double_t errData_bgr=TMath::Sqrt
// "correct" the data, using the Monte Carlo and neglecting off-diagonals
Double_t fCorr=(histBbbSignalGen->GetBinContent(i)/
Double_t data_bbb= data_bgr *fCorr;
// stat only error
Double_t errData_stat_bbb = errData*fCorr;
// stat plus background subtraction
Double_t errData_statbgr_bbb = errData_bgr*fCorr;
// estimate systematic error by repeating the exercise
// using the MC with systematic shifts applied
Double_t fCorr_sys=(histBbbSignalGenSys->GetBinContent(i)/
Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
// add systematic shift quadratically and get total error
Double_t errData_total_bbb=
// get results from real unfolding
Double_t data_unfold= histUnfoldStat->GetBinContent(i);
Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i);
Double_t errData_statbgr_unfold=histUnfoldOutput->GetBinError(i);
Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);
// compare
("%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
(" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:92
bool Bool_t
Definition RtypesCore.h:63
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kCyan
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
double pow(double, double)
R__EXTERN TStyle * gStyle
Definition TStyle.h:412
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:23
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:769
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:8903
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
Return maximum value smaller than maxval of bins in the range, unless the value has been overridden b...
Definition TH1.cxx:8390
virtual void SetMaximum(Double_t maximum=-1111)
Definition TH1.h:398
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:9046
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:399
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
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9062
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4993
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
Int_t Fill(Double_t)
Invalid Fill method.
Definition TH2.cxx:294
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:2507
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1177
Random number generator class based on M.
Definition TRandom3.h:27
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma.
Definition TRandom.cxx:225
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552
Base class for spline implementation containing the Draw/Paint methods.
Definition TSpline.h:31
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition TSpline.cxx:101
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2331
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition TStyle.cxx:1541
An algorithm to unfold distributions from detector to truth level.
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeBinWidth
scale factors from multidimensional bin width
@ kSysErrModeMatrix
matrix is an alternative to the default matrix, the errors are the difference to the original matrix
Definition TUnfoldSys.h:104
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition TUnfold.cxx:3681
type of extra constraint
Definition TUnfold.h:109
@ kEConstraintArea
enforce preservation of the area
Definition TUnfold.h:115
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
const Double_t sigma
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Double_t Exp(Double_t x)
Definition TMath.h:727
Double_t Log(Double_t x)
Definition TMath.h:760
Double_t Sqrt(Double_t x)
Definition TMath.h:691
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition TMath.h:735
static void output(int code)
Definition gifencode.c:226

Version 17.6, in parallel to changes in TUnfold


  • 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, change to use the class TUnfoldDensity
  • Version 16.1, parallel to changes in TUnfold
  • Version 16.0, parallel to changes in TUnfold
  • Version 15, simple example including background subtraction

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/.

Stefan Schmitt DESY, 14.10.2008

Definition in file testUnfold3.C.