TSVDUnfoldExample

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)

Author: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, August 17, 2022 at 09:34 AM.

Definition of a helper function:

In [ ]:
%%cpp -d

#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;
   }
}
In [ ]:
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);

TRandom3 R;

const Double_t cutdummy= -99999.0;

Data/MC toy generation

The MC input

In [ ]:
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

In [ ]:
TH1D *data = new TH1D("data", "data", nbins, -10.0, 10.0);

Data "truth" distribution to test the unfolding

In [ ]:
TH1D *datatrue = new TH1D("datatrue", "data truth", nbins, -10.0, 10.0);

Statistical covariance matrix

In [ ]:
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.

In [ ]:
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.

In [ ]:
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

In [ ]:
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

In [ ]:
TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );

It is possible to normalise unfolded spectrum to unit area

In [ ]:
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
In [ ]:
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
In [ ]:
TH1D* ddist = tsvdunf->GetD();

Get the distribution of the singular values

In [ ]:
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.

In [ ]:
TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );

Now compute the error matrix on the unfolded distribution originating from the finite detector matrix statistics

In [ ]:
TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );

Sum up the two (they are uncorrelated)

In [ ]:
ustatcov->Add( uadetcov );

Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics.

In [ ]:
TH2D* utaucov = tsvdunf->GetXtau();
utaucov->Add( uadetcov );

Get the computed inverse of the covariance matrix

In [ ]:
TH2D* uinvcov = tsvdunf->GetXinv();

Only plotting stuff below

In [ ]:
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

In [ ]:
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

In [ ]:
unfres->Draw("same");
datatrue->Draw("same");
data->Draw("same");
xini->Draw("same");

leg->Draw();

covariance matrix

In [ ]:
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

In [ ]:
TVirtualPad * c3 = c12->cd(2);
c3->SetLogy();

TLine *line = new TLine( 0.,1.,40.,1. );
line->SetLineStyle(2);

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();

Draw all canvases

In [ ]:
gROOT->GetListOfCanvases()->Draw()