{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1c8039e6",
   "metadata": {},
   "source": [
    "# TSVDUnfoldExample\n",
    " Data unfolding using Singular Value Decomposition\n",
    "\n",
    "TSVDUnfold example\n",
    "\n",
    "Data unfolding using Singular Value Decomposition (hep-ph/9509307)\n",
    "\n",
    "Example distribution and smearing model from Tim Adye (RAL)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "**Author:**  Kerstin Tackmann, Andreas Hoecker, Heiko Lacker  \n",
    "<i><small>This notebook tutorial was automatically generated with <a href= \"https://github.com/root-project/root/blob/master/documentation/doxygen/converttonotebook.py\">ROOTBOOK-izer</a> from the macro found in the ROOT repository  on Tuesday, May 19, 2026 at 08:11 PM.</small></i>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3bb857e",
   "metadata": {},
   "source": [
    " Definition of a helper function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eedbe78c",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "#include <iostream>\n",
    "\n",
    "#include \"TROOT.h\"\n",
    "#include \"TSystem.h\"\n",
    "#include \"TStyle.h\"\n",
    "#include \"TRandom3.h\"\n",
    "#include \"TString.h\"\n",
    "#include \"TMath.h\"\n",
    "#include \"TH1D.h\"\n",
    "#include \"TH2D.h\"\n",
    "#include \"TLegend.h\"\n",
    "#include \"TCanvas.h\"\n",
    "#include \"TColor.h\"\n",
    "#include \"TLine.h\"\n",
    "\n",
    "#include \"TSVDUnfold.h\"\n",
    "\n",
    "\n",
    "Double_t Reconstruct( Double_t xt, TRandom3& R )\n",
    "{\n",
    "   // apply some Gaussian smearing + bias and efficiency corrections to fake reconstruction\n",
    "   const Double_t cutdummy = -99999.0;\n",
    "   Double_t xeff = 0.3 + (1.0 - 0.3)/20.0*(xt + 10.0);  // efficiency\n",
    "   Double_t x    = R.Rndm();\n",
    "   if (x > xeff) return cutdummy;\n",
    "   else {\n",
    "     Double_t xsmear= R.Gaus(-2.5,0.2); // bias and smear\n",
    "     return xt+xsmear;\n",
    "   }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "784767dd",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gROOT->SetStyle(\"Plain\");\n",
    "gStyle->SetOptStat(0);\n",
    "\n",
    "TRandom3 R;\n",
    "\n",
    "const Double_t cutdummy= -99999.0;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "711bc951",
   "metadata": {},
   "source": [
    "--------------------------------------\n",
    "Data/MC toy generation\n",
    "\n",
    "The MC input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2af229dd",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Int_t nbins = 40;\n",
    "TH1D *xini = new TH1D(\"xini\", \"MC truth\", nbins, -10.0, 10.0);\n",
    "TH1D *bini = new TH1D(\"bini\", \"MC reco\", nbins, -10.0, 10.0);\n",
    "TH2D *Adet = new TH2D(\"Adet\", \"detector response\", nbins, -10.0, 10.0, nbins, -10.0, 10.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1004d89a",
   "metadata": {},
   "source": [
    "Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aacb0f7",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH1D *data = new TH1D(\"data\", \"data\", nbins, -10.0, 10.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7f120b1",
   "metadata": {},
   "source": [
    "Data \"truth\" distribution to test the unfolding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "894c92b8",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH1D *datatrue = new TH1D(\"datatrue\", \"data truth\", nbins, -10.0, 10.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d85dc846",
   "metadata": {},
   "source": [
    "Statistical covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "450cea36",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH2D *statcov = new TH2D(\"statcov\", \"covariance matrix\", nbins, -10.0, 10.0, nbins, -10.0, 10.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3525394",
   "metadata": {},
   "source": [
    "Fill the MC using a Breit-Wigner, mean 0.3 and width 2.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f8ba77",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (Int_t i= 0; i<100000; i++) {\n",
    "   Double_t xt = R.BreitWigner(0.3, 2.5);\n",
    "   xini->Fill(xt);\n",
    "   Double_t x = Reconstruct( xt, R );\n",
    "   if (x != cutdummy) {\n",
    "      Adet->Fill(x, xt);\n",
    "      bini->Fill(x);\n",
    "   }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e11cd4da",
   "metadata": {},
   "source": [
    "Fill the \"data\" with a Gaussian, mean 0 and width 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af6658a0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (Int_t i=0; i<10000; i++) {\n",
    "   Double_t xt = R.Gaus(0.0, 2.0);\n",
    "   datatrue->Fill(xt);\n",
    "   Double_t x = Reconstruct( xt, R );\n",
    "   if (x != cutdummy)\n",
    "   data->Fill(x);\n",
    "}\n",
    "\n",
    "cout << \"Created toy distributions and errors for: \" << endl;\n",
    "cout << \"... \\\"true MC\\\"   and \\\"reconstructed (smeared) MC\\\"\" << endl;\n",
    "cout << \"... \\\"true data\\\" and \\\"reconstructed (smeared) data\\\"\" << endl;\n",
    "cout << \"... the \\\"detector response matrix\\\"\" << endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba13e2c0",
   "metadata": {},
   "source": [
    "Fill the data covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c8a69b2",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (int i=1; i<=data->GetNbinsX(); i++) {\n",
    "    statcov->SetBinContent(i,i,data->GetBinError(i)*data->GetBinError(i));\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "862ad77a",
   "metadata": {},
   "source": [
    "----------------------------\n",
    "Here starts the actual unfolding\n",
    "\n",
    "Create TSVDUnfold object and initialise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7106029d",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "901d6c26",
   "metadata": {},
   "source": [
    "It is possible to normalise unfolded spectrum to unit area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c776748e",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tsvdunf->SetNormalize( kFALSE ); // no normalisation here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6de72893",
   "metadata": {},
   "source": [
    "Perform the unfolding with regularisation parameter kreg = 13\n",
    "- the larger kreg, the finer grained the unfolding, but the more fluctuations occur\n",
    "- the smaller kreg, the stronger is the regularisation and the bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "484fad00",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH1D* unfres = tsvdunf->Unfold( 13 );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb3556a2",
   "metadata": {},
   "source": [
    "Get the distribution of the d to cross check the regularization\n",
    "- choose kreg to be the point where |d_i| stop being statistically significantly >>1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9e809ee",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH1D* ddist = tsvdunf->GetD();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24a51ea5",
   "metadata": {},
   "source": [
    "Get the distribution of the singular values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df2a4a2d",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH1D* svdist = tsvdunf->GetSV();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "861aade5",
   "metadata": {},
   "source": [
    "Compute the error matrix for the unfolded spectrum using toy MC\n",
    "using the measured covariance matrix as input to generate the toys\n",
    "100 toys should usually be enough\n",
    "The same method can be used for different covariance matrices separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6958a2e",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1d16105",
   "metadata": {},
   "source": [
    "Now compute the error matrix on the unfolded distribution originating\n",
    "from the finite detector matrix statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e00db9",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe2b1e75",
   "metadata": {},
   "source": [
    "Sum up the two (they are uncorrelated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a23936d",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ustatcov->Add( uadetcov );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895c8a04",
   "metadata": {},
   "source": [
    "Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2f49abb",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH2D* utaucov = tsvdunf->GetXtau();\n",
    "utaucov->Add( uadetcov );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0074912a",
   "metadata": {},
   "source": [
    "Get the computed inverse of the covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f7fb37",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TH2D* uinvcov = tsvdunf->GetXinv();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "534e5c2d",
   "metadata": {},
   "source": [
    "---------------------------------\n",
    "Only plotting stuff below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4963a1e",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (int i=1; i<=unfres->GetNbinsX(); i++) {\n",
    "   unfres->SetBinError(i, TMath::Sqrt(utaucov->GetBinContent(i,i)));\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b073584f",
   "metadata": {},
   "source": [
    "Renormalize just to be able to plot on the same scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d44e5efe",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xini->Scale(0.7*datatrue->Integral()/xini->Integral());\n",
    "\n",
    "TLegend *leg = new TLegend(0.58,0.60,0.99,0.88);\n",
    "leg->SetBorderSize(0);\n",
    "leg->SetFillColor(0);\n",
    "leg->SetFillStyle(0);\n",
    "leg->AddEntry(unfres,\"Unfolded Data\",\"p\");\n",
    "leg->AddEntry(datatrue,\"True Data\",\"l\");\n",
    "leg->AddEntry(data,\"Reconstructed Data\",\"l\");\n",
    "leg->AddEntry(xini,\"True MC\",\"l\");\n",
    "\n",
    "TCanvas *c1 = new TCanvas( \"c1\", \"Unfolding toy example with TSVDUnfold\", 1000, 900 );\n",
    "\n",
    "c1->Divide(1,2);\n",
    "TVirtualPad * c11 = c1->cd(1);\n",
    "\n",
    "TH1D* frame = new TH1D( *unfres );\n",
    "frame->SetTitle( \"Unfolding toy example with TSVDUnfold\" );\n",
    "frame->GetXaxis()->SetTitle( \"x variable\" );\n",
    "frame->GetYaxis()->SetTitle( \"Events\" );\n",
    "frame->GetXaxis()->SetTitleOffset( 1.25 );\n",
    "frame->GetYaxis()->SetTitleOffset( 1.29 );\n",
    "frame->Draw();\n",
    "\n",
    "data->SetLineStyle(kDashed);\n",
    "data->SetLineColor(4);\n",
    "data->SetLineWidth(2);\n",
    "unfres->SetMarkerStyle(20);\n",
    "datatrue->SetLineColor(2);\n",
    "datatrue->SetLineWidth(2);\n",
    "xini->SetLineStyle(kDashed);\n",
    "xini->SetLineColor(8);\n",
    "xini->SetLineWidth(2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a8ef3a0",
   "metadata": {},
   "source": [
    "------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51ca82c",
   "metadata": {},
   "source": [
    "add histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee677d00",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unfres->Draw(\"same\");\n",
    "datatrue->Draw(\"same\");\n",
    "data->Draw(\"same\");\n",
    "xini->Draw(\"same\");\n",
    "\n",
    "leg->Draw();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81c0187c",
   "metadata": {},
   "source": [
    "covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86340328",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TVirtualPad * c12 = c1->cd(2);\n",
    "c12->Divide(2,1);\n",
    "TVirtualPad * c2 = c12->cd(1);\n",
    "c2->SetRightMargin   ( 0.15         );\n",
    "\n",
    "TH2D* covframe = new TH2D( *ustatcov );\n",
    "covframe->SetTitle( \"TSVDUnfold covariance matrix\" );\n",
    "covframe->GetXaxis()->SetTitle( \"x variable\" );\n",
    "covframe->GetYaxis()->SetTitle( \"x variable\" );\n",
    "covframe->GetXaxis()->SetTitleOffset( 1.25 );\n",
    "covframe->GetYaxis()->SetTitleOffset( 1.29 );\n",
    "covframe->Draw();\n",
    "\n",
    "ustatcov->SetLineWidth( 2 );\n",
    "ustatcov->Draw( \"colzsame\" );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38b1d0b2",
   "metadata": {},
   "source": [
    "distribution of the d quantity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7625ae1f",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TVirtualPad * c3 = c12->cd(2);\n",
    "c3->SetLogy();\n",
    "\n",
    "TLine *line = new TLine( 0.,1.,40.,1. );\n",
    "line->SetLineStyle(kDashed);\n",
    "\n",
    "TH1D* dframe = new TH1D( *ddist );\n",
    "dframe->SetTitle( \"TSVDUnfold |d_{i}|\" );\n",
    "dframe->GetXaxis()->SetTitle( \"i\" );\n",
    "dframe->GetYaxis()->SetTitle( \"|d_{i}|\" );\n",
    "dframe->GetXaxis()->SetTitleOffset( 1.25 );\n",
    "dframe->GetYaxis()->SetTitleOffset( 1.29 );\n",
    "dframe->SetMinimum( 0.001 );\n",
    "dframe->Draw();\n",
    "\n",
    "ddist->SetLineWidth( 2 );\n",
    "ddist->Draw( \"same\" );\n",
    "line->Draw();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7182a98",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a17993f3",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gROOT->GetListOfCanvases()->Draw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ROOT C++",
   "language": "c++",
   "name": "root"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".C",
   "mimetype": " text/x-c++src",
   "name": "c++"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}