Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
testUnfold4.C File Reference
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.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
#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;
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.,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");
}
r
ROOT::R::TRInterface & r
Definition Object.C:4
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
kWarning
const Int_t kWarning
Definition TError.h:47
gErrorIgnoreLevel
R__EXTERN Int_t gErrorIgnoreLevel
Definition TError.h:129
gStyle
R__EXTERN TStyle * gStyle
Definition TStyle.h:412
TCanvas
The Canvas class.
Definition TCanvas.h:23
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
TH1D::Reset
virtual void Reset(Option_t *option="")
Reset.
Definition TH1.cxx:10126
TH1::GetBinError
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:8903
TH1::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition TH1.cxx:3892
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
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:6663
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4993
TH2D
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:292
TH2D::Reset
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition TH2.cxx:3846
TH2::Fill
Int_t Fill(Double_t)
Invalid Fill method.
Definition TH2.cxx:294
TPad::Divide
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
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
TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition TRandom.cxx:274
TRandom::Poisson
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:402
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552
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:2331
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:1541
TUnfoldSys
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition TUnfoldSys.h:55
TUnfold::GetTUnfoldVersion
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition TUnfold.cxx:3681
TUnfold::kEConstraintArea
@ kEConstraintArea
enforce preservation of the area
Definition TUnfold.h:115
TUnfold::kEConstraintNone
@ kEConstraintNone
use no extra constraint
Definition TUnfold.h:112
TUnfold::kRegModeSize
@ kRegModeSize
regularise the amplitude of the output distribution
Definition TUnfold.h:125
TUnfold::kHistMapOutputHoriz
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition TUnfold.h:142
output
static void output(int code)
Definition gifencode.c:226

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.