Logo ROOT  
Reference Guide
testUnfold4.C File Reference

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();
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
// 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
// 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));
}
}
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");
}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
const Int_t kWarning
Definition: TError.h:38
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:105
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
The Canvas class.
Definition: TCanvas.h:27
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9741
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
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:3808
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
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:6330
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH2.cxx:3816
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
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 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:263
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
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:2311
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:1542
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition: TUnfold.cxx:3680
@ kEConstraintArea
enforce preservation of the area
Definition: TUnfold.h:115
@ kEConstraintNone
use no extra constraint
Definition: TUnfold.h:112
@ kRegModeSize
regularise the amplitude of the output distribution
Definition: TUnfold.h:125
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
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.