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 Wednesday, April 17, 2024 at 11:24 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 std::cout;
TRandom *rnd=nullptr;
input_line_43:11:10: fatal error: 'TUnfoldDensity.h' file not found #include "TUnfoldDensity.h" ^~~~~~~~~~~~~~~~~~
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;
}
input_line_44:6:15: error: use of undeclared identifier 'rnd' Double_t f=rnd->Rndm(); ^
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;
}
input_line_45:10:15: error: use of undeclared identifier 'rnd' Double_t f=rnd->Rndm(); ^ input_line_45:13:9: error: use of undeclared identifier 'rnd' r=rnd->Gaus(shapeParm[1],shapeParm[2]); ^ input_line_45:15:9: error: use of undeclared identifier 'rnd' r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3]; ^
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();
input_line_57:2:2: warning: declaration without the 'auto' keyword is deprecated: function '__cling_Un1Qu322' [-Wdeprecated-declarations] rnd=new TRandom3(); ^~~ auto
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");
input_line_65:2:31: error: 'TUnfold' is not a class, namespace, or enumeration cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n"; ^ input_line_65:2:31: note: 'TUnfold' declared here input_line_65:11:18: error: use of undeclared identifier 'GenerateGenEvent' Int_t iGen=GenerateGenEvent(nGen,genFrac); ^ input_line_65:12:21: error: use of undeclared identifier 'GenerateRecEvent' Double_t yObs=GenerateRecEvent(genShape[iGen]); ^ input_line_65:18:18: error: use of undeclared identifier 'GenerateGenEvent' Int_t iGen=GenerateGenEvent(nGen,genFrac); ^ input_line_65:19:21: error: use of undeclared identifier 'GenerateRecEvent' Double_t yObs=GenerateRecEvent(genShape[iGen]); ^ input_line_65:30:4: error: unknown type name 'TUnfoldSys' TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz, ^ input_line_65:30:35: error: use of undeclared identifier 'TUnfold' TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz, ^ input_line_65:31:22: error: use of undeclared identifier 'TUnfold' TUnfold::kRegModeSize,TUnfold::kEConstraintNone); ^ input_line_65:31:44: error: use of undeclared identifier 'TUnfold' TUnfold::kRegModeSize,TUnfold::kEConstraintNone); ^ input_line_65:47:25: error: use of undeclared identifier 'TUnfold' unfold.SetConstraint(TUnfold::kEConstraintArea); ^ input_line_65:11:35: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Int_t *' (aka 'const int *') Int_t iGen=GenerateGenEvent(nGen,genFrac); ^~~~ input_line_65:11:40: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Double_t (*)[3]' (aka 'const double (*)[3]') Int_t iGen=GenerateGenEvent(nGen,genFrac); ^~~~~~~ input_line_65:12:38: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Double_t (*)[3][5]' (aka 'const double (*)[3][5]') Double_t yObs=GenerateRecEvent(genShape[iGen]); ^~~~~~~~ input_line_65:18:35: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Int_t *' (aka 'const int *') Int_t iGen=GenerateGenEvent(nGen,genFrac); ^~~~ input_line_65:18:40: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Double_t (*)[3]' (aka 'const double (*)[3]') Int_t iGen=GenerateGenEvent(nGen,genFrac); ^~~~~~~ input_line_65:19:38: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const Double_t (*)[3][5]' (aka 'const double (*)[3][5]') Double_t yObs=GenerateRecEvent(genShape[iGen]); ^~~~~~~~
Draw all canvases
gROOT->GetListOfCanvases()->Draw()
In module 'Cling_Runtime_Extra': /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<int>::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); ^