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.
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 Sunday, February 05, 2023 at 11:21 AM.
%%cpp -d
#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 namespace std;
TRandom *rnd=0;
Definition of a helper function:
%%cpp -d
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;
}
Definition of a helper function:
%%cpp -d
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;
}
Arguments are defined.
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.,0,0,0);
// 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.,0,0,0);
// 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");
TUnfold version is V17.6 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 FCN=17.3145 FROM MIGRAD STATUS=CONVERGED 59 CALLS 60 TOTAL EDM=1.33067e-10 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.53639e+02 6.05704e+00 1.02575e-02 3.01728e-06 2 Mean -1.21390e-01 3.38684e-02 6.92568e-05 -1.72950e-05 3 Sigma 1.02094e+00 2.43696e-02 1.33049e-05 2.18204e-03 FCN=8.19808 FROM MIGRAD STATUS=CONVERGED 59 CALLS 60 TOTAL EDM=6.85034e-10 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.49698e+02 5.79213e+00 7.06055e-03 3.51186e-06 2 Mean -2.59551e-01 3.46939e-02 5.10123e-05 1.03983e-04 3 Sigma 1.05666e+00 2.43573e-02 9.38061e-06 5.79231e-03 FCN=59.2962 FROM MIGRAD STATUS=CONVERGED 60 CALLS 61 TOTAL EDM=7.24941e-07 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.27750e+02 5.41128e+00 1.64364e-02 -2.21211e-04 2 Mean -4.82115e-01 4.90259e-02 1.46715e-04 -1.83701e-02 3 Sigma 1.09744e+00 3.08576e-02 2.35052e-05 -1.11293e-01 FCN=24.0114 FROM MIGRAD STATUS=CONVERGED 59 CALLS 60 TOTAL EDM=2.61633e-10 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.47374e+02 5.80099e+00 1.15432e-02 -2.73138e-06 2 Mean 1.89896e-01 3.51550e-02 8.46764e-05 2.01949e-04 3 Sigma 1.05773e+00 2.51629e-02 1.57038e-05 -3.21527e-03 FCN=17.3474 FROM MIGRAD STATUS=CONVERGED 59 CALLS 60 TOTAL EDM=3.1205e-08 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.49983e+02 5.77488e+00 1.00379e-02 -3.84601e-05 2 Mean 2.01446e-01 3.42777e-02 7.12868e-05 3.82112e-03 3 Sigma 1.04383e+00 2.32937e-02 1.28346e-05 -3.31497e-02 FCN=66.0535 FROM MIGRAD STATUS=CONVERGED 69 CALLS 70 TOTAL EDM=1.81967e-10 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 Constant 1.24001e+02 5.38113e+00 1.68582e-02 3.09845e-06 2 Mean -4.12797e-01 4.98370e-02 1.59404e-04 -2.92080e-04 3 Sigma 1.12618e+00 3.36751e-02 2.62795e-05 -5.54243e-04
Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000. Warning in <TUnfoldSys::SetInput>: One input bin has zero error, 1/error set to 1.000000.
Draw all canvases
gROOT->GetListOfCanvases()->Draw()