/// \file /// \ingroup tutorial_math /// \notebook /// 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) /// /// \macro_image /// \macro_code /// /// \authors Kerstin Tackmann, Andreas Hoecker, Heiko Lacker #include #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"); gStyle->SetOptStat(0); TRandom3 R; 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. ); 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(); }