Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
testUnfold4.C File Reference
Tutorials » Data analysis tutorials » TUnfold tutorials

Detailed Description

View in nbviewer Open in SWAN
Test program for the class TUnfoldSys.

Simple toy tests of the TUnfold package

Pseudo data (5000 events) are unfolded into three components The unfolding is performed once without and once with area constraint

Ideally, the pulls may show that the result is biased if no constraint is applied. This is expected because the true data errors are not known, and instead the sqrt(data) errors are used.

TUnfold version is V17.9
toy iteration: 0
toy iteration: 10
toy iteration: 20
toy iteration: 30
toy iteration: 40
toy iteration: 50
toy iteration: 60
toy iteration: 70
toy iteration: 80
toy iteration: 90
toy iteration: 100
toy iteration: 110
toy iteration: 120
toy iteration: 130
toy iteration: 140
toy iteration: 150
toy iteration: 160
toy iteration: 170
toy iteration: 180
toy iteration: 190
toy iteration: 200
toy iteration: 210
toy iteration: 220
toy iteration: 230
toy iteration: 240
toy iteration: 250
toy iteration: 260
toy iteration: 270
toy iteration: 280
toy iteration: 290
toy iteration: 300
toy iteration: 310
toy iteration: 320
toy iteration: 330
toy iteration: 340
toy iteration: 350
toy iteration: 360
toy iteration: 370
toy iteration: 380
toy iteration: 390
toy iteration: 400
toy iteration: 410
toy iteration: 420
toy iteration: 430
toy iteration: 440
toy iteration: 450
toy iteration: 460
toy iteration: 470
toy iteration: 480
toy iteration: 490
toy iteration: 500
toy iteration: 510
toy iteration: 520
toy iteration: 530
toy iteration: 540
toy iteration: 550
toy iteration: 560
toy iteration: 570
toy iteration: 580
toy iteration: 590
toy iteration: 600
toy iteration: 610
toy iteration: 620
toy iteration: 630
toy iteration: 640
toy iteration: 650
toy iteration: 660
toy iteration: 670
toy iteration: 680
toy iteration: 690
toy iteration: 700
toy iteration: 710
toy iteration: 720
toy iteration: 730
toy iteration: 740
toy iteration: 750
toy iteration: 760
toy iteration: 770
toy iteration: 780
toy iteration: 790
toy iteration: 800
toy iteration: 810
toy iteration: 820
toy iteration: 830
toy iteration: 840
toy iteration: 850
toy iteration: 860
toy iteration: 870
toy iteration: 880
toy iteration: 890
toy iteration: 900
toy iteration: 910
toy iteration: 920
toy iteration: 930
toy iteration: 940
toy iteration: 950
toy iteration: 960
toy iteration: 970
toy iteration: 980
toy iteration: 990
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 17.3145
NDf = 12
Edm = 7.81966e-06
NCalls = 52
Constant = 153.648 +/- 6.05711
Mean = -0.121307 +/- 0.0338681
Sigma = 1.0209 +/- 0.0243688 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 8.19808
NDf = 12
Edm = 2.68037e-06
NCalls = 53
Constant = 149.692 +/- 5.79204
Mean = -0.259587 +/- 0.0346952
Sigma = 1.05667 +/- 0.0243581 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 59.2962
NDf = 12
Edm = 3.6247e-07
NCalls = 61
Constant = 127.75 +/- 5.41128
Mean = -0.482115 +/- 0.0490259
Sigma = 1.09744 +/- 0.0308576 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 24.0114
NDf = 12
Edm = 3.60094e-06
NCalls = 53
Constant = 147.382 +/- 5.80116
Mean = 0.189905 +/- 0.0351528
Sigma = 1.05774 +/- 0.0251621 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 17.3474
NDf = 12
Edm = 7.55094e-06
NCalls = 53
Constant = 149.977 +/- 5.77475
Mean = 0.201466 +/- 0.034281
Sigma = 1.04379 +/- 0.0232952 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 66.0535
NDf = 12
Edm = 1.27624e-06
NCalls = 61
Constant = 124.004 +/- 5.38123
Mean = -0.412836 +/- 0.0498374
Sigma = 1.12618 +/- 0.0336756 (limited)
#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 <TError.h>
#include "TUnfoldDensity.h"
using std::cout;
TRandom *rnd=nullptr;
Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) {
// choose an integer random number in the range [0,nmax]
// (the generator level bin)
// depending on the probabilities
// probability[0],probability[1],...probability[nmax-1]
Double_t f=rnd->Rndm();
Int_t r=0;
while((r<nmax)&&(f>=probability[r])) {
f -= probability[r];
r++;
}
return r;
}
Double_t GenerateRecEvent(const Double_t *shapeParm) {
// return a coordinate (the reconstructed variable)
// depending on shapeParm[]
// shapeParm[0]: fraction of events with Gaussian distribution
// shapeParm[1]: mean of Gaussian
// shapeParm[2]: width of Gaussian
// (1-shapeParm[0]): fraction of events with flat distribution
// shapeParm[3]: minimum of flat component
// shapeParm[4]: maximum of flat component
Double_t f=rnd->Rndm();
Double_t r;
if(f<shapeParm[0]) {
r=rnd->Gaus(shapeParm[1],shapeParm[2]);
} else {
r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3];
}
return r;
}
void testUnfold4(bool printInfo = false)
{
// switch off printing Info messages
if (!printInfo) gErrorIgnoreLevel = kWarning;
// switch on histogram errors
TH1::SetDefaultSumw2();
// random generator
rnd=new TRandom3();
// data and MC number of events
Double_t const nData0= 500.0;
Double_t const nMC0 = 50000.0;
// Binning
// reconstructed variable (0-10)
Int_t const nDet=15;
Double_t const xminDet=0.0;
Double_t const xmaxDet=15.0;
// signal binning (three shapes: 0,1,2)
Int_t const nGen=3;
Double_t const xminGen=-0.5;
Double_t const xmaxGen= 2.5;
// parameters
// fraction of events per signal shape
static const Double_t genFrac[]={0.3,0.6,0.1};
// signal shapes
static const Double_t genShape[][5]=
{{1.0,2.0,1.5,0.,15.},
{1.0,7.0,2.5,0.,15.},
{0.0,0.0,0.0,0.,15.}};
// define DATA histograms
// observed data distribution
TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet);
// define MC histograms
// matrix of migrations
TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)",
nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen);
TH1D **histPullNC=new TH1D* [nGen];
TH1D **histPullArea=new TH1D* [nGen];
for(int i=0;i<nGen;i++) {
histPullNC[i]=new TH1D(TString::Format("PullNC%d",i),"pull",15,-3.,3.);
histPullArea[i]=new TH1D(TString::Format("PullArea%d",i),"pull",15,-3.,3.);
}
// this method is new in version 16 of TUnfold
cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
for(int itoy=0;itoy<1000;itoy++) {
if(!(itoy %10)) cout<<"toy iteration: "<<itoy<<"\n";
histDetDATA->Reset();
histGenDetMC->Reset();
Int_t nData=rnd->Poisson(nData0);
for(Int_t i=0;i<nData;i++) {
Int_t iGen=GenerateGenEvent(nGen,genFrac);
Double_t yObs=GenerateRecEvent(genShape[iGen]);
histDetDATA->Fill(yObs);
}
Int_t nMC=rnd->Poisson(nMC0);
for(Int_t i=0;i<nMC;i++) {
Int_t iGen=GenerateGenEvent(nGen,genFrac);
Double_t yObs=GenerateRecEvent(genShape[iGen]);
histGenDetMC->Fill(iGen,yObs);
}
/* for(Int_t ix=0;ix<=histGenDetMC->GetNbinsX()+1;ix++) {
for(Int_t iy=0;iy<=histGenDetMC->GetNbinsY()+1;iy++) {
cout<<ix<<iy<<" : "<<histGenDetMC->GetBinContent(ix,iy)<<"\n";
}
} */
//========================
// unfolding
TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz,
TUnfold::kRegModeSize,TUnfold::kEConstraintNone);
// define the input vector (the measured data distribution)
unfold.SetInput(histDetDATA,0.0,1.0);
// run the unfolding
unfold.ScanLcurve(50,0.,0.,nullptr,nullptr,nullptr);
// fill pull distributions without constraint
unfold.GetOutput(histUnfold);
for(int i=0;i<nGen;i++) {
histPullNC[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
histUnfold->GetBinError(i+1));
}
// repeat unfolding on the same data, now with Area constraint
unfold.SetConstraint(TUnfold::kEConstraintArea);
// run the unfolding
unfold.ScanLcurve(50,0.,0.,nullptr,nullptr,nullptr);
// fill pull distributions with constraint
unfold.GetOutput(histUnfold);
for(int i=0;i<nGen;i++) {
histPullArea[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
histUnfold->GetBinError(i+1));
}
}
TCanvas output;
output.Divide(3,2);
gStyle->SetOptFit(1111);
for(int i=0;i<nGen;i++) {
output.cd(i+1);
histPullNC[i]->Fit("gaus");
histPullNC[i]->Draw();
}
for(int i=0;i<nGen;i++) {
output.cd(i+4);
histPullArea[i]->Fit("gaus");
histPullArea[i]->Draw();
}
output.SaveAs("testUnfold4.ps");
}
f
#define f(i)
Definition RSha256.hxx:104
Int_t
int Int_t
Definition RtypesCore.h:45
Double_t
double Double_t
Definition RtypesCore.h:59
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:358
kWarning
constexpr Int_t kWarning
Definition TError.h:46
gErrorIgnoreLevel
Int_t gErrorIgnoreLevel
Error handling routines.
Definition TError.cxx:31
r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Definition TGWin32VirtualXProxy.cxx:168
gStyle
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
TCanvas
The Canvas class.
Definition TCanvas.h:23
TH1D
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:693
TH1::SetDefaultSumw2
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:6701
TH2D
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:351
TRandom3
Random number generator class based on M.
Definition TRandom3.h:27
TRandom
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
TString::Format
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:2378
TStyle::SetOptFit
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:1593
TUnfoldSys
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition TUnfoldSys.h:59
TUnfold::GetTUnfoldVersion
static const char * GetTUnfoldVersion(void)
return a string describing the TUnfold version
Definition TUnfold.cxx:3717
TUnfold::kEConstraintArea
@ kEConstraintArea
enforce preservation of the area
Definition TUnfold.h:119
TUnfold::kEConstraintNone
@ kEConstraintNone
use no extra constraint
Definition TUnfold.h:116
TUnfold::kRegModeSize
@ kRegModeSize
regularise the amplitude of the output distribution
Definition TUnfold.h:129
TUnfold::kHistMapOutputHoriz
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition TUnfold.h:146
output
static void output()

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 16.1, parallel to changes in TUnfold
  • Version 16.0, parallel to changes in TUnfold
  • Version 15, use L-curve scan to scan the average correlation

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