Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TSVDUnfoldExample.C File Reference

Detailed Description

View in nbviewer Open in SWAN
Data unfolding using Singular Value Decomposition

TSVDUnfold example

Data unfolding using Singular Value Decomposition (hep-ph/9509307)

Example distribution and smearing model from Tim Adye (RAL)

#include <iostream>
#include "TROOT.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TString.h"
#include "TMath.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TLine.h"
#include "TSVDUnfold.h"
Double_t Reconstruct( Double_t xt, TRandom3& R )
{
// apply some Gaussian smearing + bias and efficiency corrections to fake reconstruction
const Double_t cutdummy = -99999.0;
Double_t xeff = 0.3 + (1.0 - 0.3)/20.0*(xt + 10.0); // efficiency
Double_t x = R.Rndm();
if (x > xeff) return cutdummy;
else {
Double_t xsmear= R.Gaus(-2.5,0.2); // bias and smear
return xt+xsmear;
}
}
void TSVDUnfoldExample()
{
gROOT->SetStyle("Plain");
const Double_t cutdummy= -99999.0;
// --------------------------------------
// Data/MC toy generation
//
// The MC input
Int_t nbins = 40;
TH1D *xini = new TH1D("xini", "MC truth", nbins, -10.0, 10.0);
TH1D *bini = new TH1D("bini", "MC reco", nbins, -10.0, 10.0);
TH2D *Adet = new TH2D("Adet", "detector response", nbins, -10.0, 10.0, nbins, -10.0, 10.0);
// Data
TH1D *data = new TH1D("data", "data", nbins, -10.0, 10.0);
// Data "truth" distribution to test the unfolding
TH1D *datatrue = new TH1D("datatrue", "data truth", nbins, -10.0, 10.0);
// Statistical covariance matrix
TH2D *statcov = new TH2D("statcov", "covariance matrix", nbins, -10.0, 10.0, nbins, -10.0, 10.0);
// Fill the MC using a Breit-Wigner, mean 0.3 and width 2.5.
for (Int_t i= 0; i<100000; i++) {
Double_t xt = R.BreitWigner(0.3, 2.5);
xini->Fill(xt);
Double_t x = Reconstruct( xt, R );
if (x != cutdummy) {
Adet->Fill(x, xt);
bini->Fill(x);
}
}
// Fill the "data" with a Gaussian, mean 0 and width 2.
for (Int_t i=0; i<10000; i++) {
Double_t xt = R.Gaus(0.0, 2.0);
datatrue->Fill(xt);
Double_t x = Reconstruct( xt, R );
if (x != cutdummy)
data->Fill(x);
}
cout << "Created toy distributions and errors for: " << endl;
cout << "... \"true MC\" and \"reconstructed (smeared) MC\"" << endl;
cout << "... \"true data\" and \"reconstructed (smeared) data\"" << endl;
cout << "... the \"detector response matrix\"" << endl;
// Fill the data covariance matrix
for (int i=1; i<=data->GetNbinsX(); i++) {
statcov->SetBinContent(i,i,data->GetBinError(i)*data->GetBinError(i));
}
// ----------------------------
// Here starts the actual unfolding
//
// Create TSVDUnfold object and initialise
TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );
// It is possible to normalise unfolded spectrum to unit area
tsvdunf->SetNormalize( kFALSE ); // no normalisation here
// Perform the unfolding with regularisation parameter kreg = 13
// - the larger kreg, the finer grained the unfolding, but the more fluctuations occur
// - the smaller kreg, the stronger is the regularisation and the bias
TH1D* unfres = tsvdunf->Unfold( 13 );
// Get the distribution of the d to cross check the regularization
// - choose kreg to be the point where |d_i| stop being statistically significantly >>1
TH1D* ddist = tsvdunf->GetD();
// Get the distribution of the singular values
TH1D* svdist = tsvdunf->GetSV();
// Compute the error matrix for the unfolded spectrum using toy MC
// using the measured covariance matrix as input to generate the toys
// 100 toys should usually be enough
// The same method can be used for different covariance matrices separately.
TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );
// Now compute the error matrix on the unfolded distribution originating
// from the finite detector matrix statistics
TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );
// Sum up the two (they are uncorrelated)
ustatcov->Add( uadetcov );
//Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics.
TH2D* utaucov = tsvdunf->GetXtau();
utaucov->Add( uadetcov );
//Get the computed inverse of the covariance matrix
TH2D* uinvcov = tsvdunf->GetXinv();
// ---------------------------------
// Only plotting stuff below
for (int i=1; i<=unfres->GetNbinsX(); i++) {
unfres->SetBinError(i, TMath::Sqrt(utaucov->GetBinContent(i,i)));
}
// Renormalize just to be able to plot on the same scale
xini->Scale(0.7*datatrue->Integral()/xini->Integral());
TLegend *leg = new TLegend(0.58,0.60,0.99,0.88);
leg->SetBorderSize(0);
leg->SetFillColor(0);
leg->SetFillStyle(0);
leg->AddEntry(unfres,"Unfolded Data","p");
leg->AddEntry(datatrue,"True Data","l");
leg->AddEntry(data,"Reconstructed Data","l");
leg->AddEntry(xini,"True MC","l");
TCanvas *c1 = new TCanvas( "c1", "Unfolding toy example with TSVDUnfold", 1000, 900 );
c1->Divide(1,2);
TVirtualPad * c11 = c1->cd(1);
TH1D* frame = new TH1D( *unfres );
frame->SetTitle( "Unfolding toy example with TSVDUnfold" );
frame->GetXaxis()->SetTitle( "x variable" );
frame->GetYaxis()->SetTitle( "Events" );
frame->GetXaxis()->SetTitleOffset( 1.25 );
frame->GetYaxis()->SetTitleOffset( 1.29 );
frame->Draw();
data->SetLineStyle(2);
data->SetLineColor(4);
data->SetLineWidth(2);
unfres->SetMarkerStyle(20);
datatrue->SetLineColor(2);
datatrue->SetLineWidth(2);
xini->SetLineStyle(2);
xini->SetLineColor(8);
xini->SetLineWidth(2);
// ------------------------------------------------------------
// add histograms
unfres->Draw("same");
datatrue->Draw("same");
data->Draw("same");
xini->Draw("same");
leg->Draw();
// covariance matrix
TVirtualPad * c12 = c1->cd(2);
c12->Divide(2,1);
TVirtualPad * c2 = c12->cd(1);
c2->SetRightMargin ( 0.15 );
TH2D* covframe = new TH2D( *ustatcov );
covframe->SetTitle( "TSVDUnfold covariance matrix" );
covframe->GetXaxis()->SetTitle( "x variable" );
covframe->GetYaxis()->SetTitle( "x variable" );
covframe->GetXaxis()->SetTitleOffset( 1.25 );
covframe->GetYaxis()->SetTitleOffset( 1.29 );
covframe->Draw();
ustatcov->SetLineWidth( 2 );
ustatcov->Draw( "colzsame" );
// distribution of the d quantity
TVirtualPad * c3 = c12->cd(2);
c3->SetLogy();
TLine *line = new TLine( 0.,1.,40.,1. );
TH1D* dframe = new TH1D( *ddist );
dframe->SetTitle( "TSVDUnfold |d_{i}|" );
dframe->GetXaxis()->SetTitle( "i" );
dframe->GetYaxis()->SetTitle( "|d_{i}|" );
dframe->GetXaxis()->SetTitleOffset( 1.25 );
dframe->GetYaxis()->SetTitleOffset( 1.29 );
dframe->SetMinimum( 0.001 );
dframe->Draw();
ddist->SetLineWidth( 2 );
ddist->Draw( "same" );
line->Draw();
}
#define R(a, b, c, d, e, f, g, h, i)
Definition RSha256.hxx:110
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
double Double_t
Definition RtypesCore.h:59
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
#define gROOT
Definition TROOT.h:406
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
virtual void SetRightMargin(Float_t rightmargin)
Set Pad right margin in fraction of the pad width.
Definition TAttPad.cxx:119
The Canvas class.
Definition TCanvas.h:23
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:670
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6739
TAxis * GetXaxis()
Definition TH1.h:324
virtual Int_t GetNbinsX() const
Definition TH1.h:297
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:826
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition TH1.cxx:9227
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3344
TAxis * GetYaxis()
Definition TH1.h:325
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:405
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition TH1.cxx:7962
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6625
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:357
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2616
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:93
Int_t Fill(Double_t) override
Invalid Fill method.
Definition TH2.cxx:393
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
Use the TLine constructor to create a simple line.
Definition TLine.h:22
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:292
Random number generator class based on M.
Definition TRandom3.h:27
SVD Approach to Data Unfolding.
Definition TSVDUnfold.h:46
TH1D * GetSV() const
Returns singular values vector.
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
void SetNormalize(Bool_t normalize)
Definition TSVDUnfold.h:66
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
TH2D * GetUnfoldCovMatrix(const TH2D *cov, Int_t ntoys, Int_t seed=1)
Determine for given input error matrix covariance matrix of unfolded spectrum from toy simulation giv...
TH2D * GetAdetCovMatrix(Int_t ntoys, Int_t seed=1)
Determine covariance matrix of unfolded spectrum from finite statistics in response matrix using pseu...
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1640
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition TVirtualPad.h:51
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)=0
virtual void SetLogy(Int_t value=1)=0
TLine * line
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
leg
Definition legend1.C:34
return c2
Definition legend2.C:14
return c3
Definition legend3.C:15
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Authors
Kerstin Tackmann, Andreas Hoecker, Heiko Lacker

Definition in file TSVDUnfoldExample.C.