{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c37fa810",
   "metadata": {},
   "source": [
    "# testUnfold5d\n",
    "Test program for the classes TUnfoldDensity and TUnfoldBinning.\n",
    "\n",
    "A toy test of the TUnfold package\n",
    "\n",
    "This is an example of unfolding a two-dimensional distribution\n",
    "also using an auxiliary measurement to constrain some background\n",
    "\n",
    "The example comprises several macros\n",
    "  testUnfold5a.C   create root files with TTree objects for\n",
    "                     signal, background and data\n",
    "           -> write files  testUnfold5_signal.root\n",
    "                           testUnfold5_background.root\n",
    "                           testUnfold5_data.root\n",
    "\n",
    "  testUnfold5b.C   create a root file with the TUnfoldBinning objects\n",
    "           -> write file  testUnfold5_binning.root\n",
    "\n",
    "  testUnfold5c.C   loop over trees and fill histograms based on the\n",
    "                     TUnfoldBinning objects\n",
    "           -> read  testUnfold5_binning.root\n",
    "                    testUnfold5_signal.root\n",
    "                    testUnfold5_background.root\n",
    "                    testUnfold5_data.root\n",
    "\n",
    "           -> write testUnfold5_histograms.root\n",
    "\n",
    "  testUnfold5d.C   run the unfolding\n",
    "           -> read  testUnfold5_histograms.root\n",
    "           -> write testUnfold5_result.root\n",
    "                    testUnfold5_result.ps\n",
    "\n",
    "\n",
    " **Version 17.6, in parallel to changes in TUnfold**\n",
    "\n",
    "#### History:\n",
    " - Version 17.5, in parallel to changes in TUnfold\n",
    " - Version 17.4, in parallel to changes in TUnfold\n",
    " - Version 17.3, in parallel to changes in TUnfold\n",
    " - Version 17.2, in parallel to changes in TUnfold\n",
    " - Version 17.1, in parallel to changes in TUnfold\n",
    " - Version 17.0 example for multi-dimensional unfolding\n",
    "\n",
    " This file is part of TUnfold.\n",
    "\n",
    " TUnfold is free software: you can redistribute it and/or modify\n",
    " it under the terms of the GNU General Public License as published by\n",
    " the Free Software Foundation, either version 3 of the License, or\n",
    " (at your option) any later version.\n",
    "\n",
    " TUnfold is distributed in the hope that it will be useful,\n",
    " but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n",
    " GNU General Public License for more details.\n",
    "\n",
    " You should have received a copy of the GNU General Public License\n",
    " along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.\n",
    "\n",
    "\n",
    "\n",
    "**Author:** Stefan Schmitt DESY, 14.10.2008  \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": "code",
   "execution_count": null,
   "id": "63a0ea13",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "#include <iostream>\n",
    "#include <cmath>\n",
    "#include <TMath.h>\n",
    "#include <TCanvas.h>\n",
    "#include <TStyle.h>\n",
    "#include <TGraph.h>\n",
    "#include <TFile.h>\n",
    "#include <TH1.h>\n",
    "#include \"TUnfoldDensity.h\"\n",
    "\n",
    "using std::cout;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75926cb7",
   "metadata": {},
   "source": [
    "switch on histogram errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5149873a",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH1::SetDefaultSumw2();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0424db9b",
   "metadata": {},
   "source": [
    "==============================================\n",
    "step 1 : open output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56d8b803",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TFile *outputFile=new TFile(\"testUnfold5_results.root\",\"recreate\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1282eb3",
   "metadata": {},
   "source": [
    "==============================================\n",
    "step 2 : read binning schemes and input histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2daccf84",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TFile *inputFile=new TFile(\"testUnfold5_histograms.root\");\n",
    "\n",
    " outputFile->cd();\n",
    "\n",
    " TUnfoldBinning *detectorBinning,*generatorBinning;\n",
    "\n",
    " inputFile->GetObject(\"detector\",detectorBinning);\n",
    " inputFile->GetObject(\"generator\",generatorBinning);\n",
    "\n",
    " if((!detectorBinning)||(!generatorBinning)) {\n",
    "    cout<<\"problem to read binning schemes\\n\";\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748f9015",
   "metadata": {},
   "source": [
    "save binning schemes to output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14a65ede",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " detectorBinning->Write();\n",
    " generatorBinning->Write();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5241c9f",
   "metadata": {},
   "source": [
    "read histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd47f838",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH1 *histDataReco,*histDataTruth;\n",
    " TH2 *histMCGenRec;\n",
    "\n",
    " inputFile->GetObject(\"histDataReco\",histDataReco);\n",
    " inputFile->GetObject(\"histDataTruth\",histDataTruth);\n",
    " inputFile->GetObject(\"histMCGenRec\",histMCGenRec);\n",
    "\n",
    "#ifdef TEST_ZERO_UNCORR_ERROR"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bdb24b8",
   "metadata": {},
   "source": [
    "special test (bug in version 17.2 and below)\n",
    "set all errors in hisMCGenRec to zero\n",
    "-> program will crash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eb0fa87",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " for(int i=0;i<=histMCGenRec->GetNbinsX()+1;i++) {\n",
    "    for(int j=0;j<=histMCGenRec->GetNbinsY()+1;j++) {\n",
    "       histMCGenRec->SetBinError(i,j,0.0);\n",
    "    }\n",
    " }\n",
    "#endif\n",
    "\n",
    " histDataReco->Write();\n",
    " histDataTruth->Write();\n",
    " histMCGenRec->Write();\n",
    "\n",
    " if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {\n",
    "    cout<<\"problem to read input histograms\\n\";\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f681b207",
   "metadata": {},
   "source": [
    "========================\n",
    "Step 3: unfolding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "998172f8",
   "metadata": {},
   "source": [
    "preserve the area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06f8da23",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2606cd0",
   "metadata": {},
   "source": [
    "basic choice of regularisation scheme:\n",
    "curvature (second derivative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b449da9d",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1626e85",
   "metadata": {},
   "source": [
    "density flags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d629d0d3",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TUnfoldDensity::EDensityMode densityFlags=\n",
    "   TUnfoldDensity::kDensityModeBinWidth;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0fb2604",
   "metadata": {},
   "source": [
    "detailed steering for regularisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f9aebb9",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " const char *REGULARISATION_DISTRIBUTION=nullptr;\n",
    " const char *REGULARISATION_AXISSTEERING=\"*[B]\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efa162e7",
   "metadata": {},
   "source": [
    "set up matrix of migrations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77101831",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz,\n",
    "                       regMode,constraintMode,densityFlags,\n",
    "                       generatorBinning,detectorBinning,\n",
    "           REGULARISATION_DISTRIBUTION,\n",
    "           REGULARISATION_AXISSTEERING);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f49600",
   "metadata": {},
   "source": [
    "define the input vector (the measured data distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f9dcffd",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ifdef TEST_INPUT_COVARIANCE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "229bbf70",
   "metadata": {},
   "source": [
    "special test to use input covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd5913bb",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH2D *inputEmatrix=\n",
    "    detectorBinning->CreateErrorMatrixHistogram(\"input_covar\",true);\n",
    " for(int i=1;i<=inputEmatrix->GetNbinsX();i++) {\n",
    "    Double_t e=histDataReco->GetBinError(i);\n",
    "    inputEmatrix->SetBinContent(i,i,e*e);\n",
    "    // test: non-zero covariance where variance is zero\n",
    "    //   if(e<=0.) inputEmatrix->SetBinContent(i,i+1,1.0);\n",
    " }\n",
    " unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);\n",
    "#else\n",
    " unfold.SetInput(histDataReco /* ,0.0,1.0 */);\n",
    "#endif"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60aadbdf",
   "metadata": {},
   "source": [
    "print matrix of regularisation conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be97e30b",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ifdef PRINT_MATRIX_L\n",
    " TH2 *histL= unfold.GetL(\"L\");\n",
    " for(Int_t j=1;j<=histL->GetNbinsY();j++) {\n",
    "    cout<<\"L[\"<<unfold.GetLBinning()->GetBinName(j)<<\"]\";\n",
    "    for(Int_t i=1;i<=histL->GetNbinsX();i++) {\n",
    "       Double_t c=histL->GetBinContent(i,j);\n",
    "       if(c!=0.0) cout<<\" [\"<<i<<\"]=\"<<c;\n",
    "    }\n",
    "    cout<<\"\\n\";\n",
    " }\n",
    "#endif"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec02bbcf",
   "metadata": {},
   "source": [
    "run the unfolding\n",
    "\n",
    "here, tau is determined by scanning the global correlation coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "393e1b45",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " Int_t nScan=30;\n",
    " TSpline *rhoLogTau=nullptr;\n",
    " TGraph *lCurve=nullptr;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8ca667d",
   "metadata": {},
   "source": [
    "for determining tau, scan the correlation coefficients\n",
    "correlation coefficients may be probed for all distributions\n",
    "or only for selected distributions\n",
    "underflow/overflow bins may be included/excluded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933c59ab",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " const char *SCAN_DISTRIBUTION=\"signal\";\n",
    " const char *SCAN_AXISSTEERING=nullptr;\n",
    "\n",
    " Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,\n",
    "                            TUnfoldDensity::kEScanTauRhoMax,\n",
    "                            SCAN_DISTRIBUTION,SCAN_AXISSTEERING,\n",
    "                            &lCurve);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24553eef",
   "metadata": {},
   "source": [
    "create graphs with one point to visualize best choice of tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01f178e1",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " Double_t t[1],rho[1],x[1],y[1];\n",
    " rhoLogTau->GetKnot(iBest,t[0],rho[0]);\n",
    " lCurve->GetPoint(iBest,x[0],y[0]);\n",
    " TGraph *bestRhoLogTau=new TGraph(1,t,rho);\n",
    " TGraph *bestLCurve=new TGraph(1,x,y);\n",
    " Double_t *tAll=new Double_t[nScan],*rhoAll=new Double_t[nScan];\n",
    " for(Int_t i=0;i<nScan;i++) {\n",
    "    rhoLogTau->GetKnot(i,tAll[i],rhoAll[i]);\n",
    " }\n",
    " TGraph *knots=new TGraph(nScan,tAll,rhoAll);\n",
    "\n",
    " cout<<\"chi**2=\"<<unfold.GetChi2A()<<\"+\"<<unfold.GetChi2L()\n",
    "     <<\" / \"<<unfold.GetNdf()<<\"\\n\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab87c217",
   "metadata": {},
   "source": [
    "===========================\n",
    "Step 4: retrieve and plot unfolding results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e8ee5be",
   "metadata": {},
   "source": [
    "get unfolding output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eda60fe",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH1 *histDataUnfold=unfold.GetOutput(\"unfolded signal\",nullptr,nullptr,nullptr,kFALSE);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ca20f4a",
   "metadata": {},
   "source": [
    "get Monte Carlo reconstructed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abc51795",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH1 *histMCReco=histMCGenRec->ProjectionY(\"histMCReco\",0,-1,\"e\");\n",
    " TH1 *histMCTruth=histMCGenRec->ProjectionX(\"histMCTruth\",0,-1,\"e\");\n",
    " Double_t scaleFactor=histDataTruth->GetSumOfWeights()/\n",
    "    histMCTruth->GetSumOfWeights();\n",
    " histMCReco->Scale(scaleFactor);\n",
    " histMCTruth->Scale(scaleFactor);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee25f701",
   "metadata": {},
   "source": [
    "get matrix of probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c6044c7",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " TH2 *histProbability=unfold.GetProbabilityMatrix(\"histProbability\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d5fe2d2",
   "metadata": {},
   "source": [
    "get global correlation coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20dbce2f",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " /* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal(\"histGlobalCorr\",nullptr,nullptr,nullptr,kFALSE);\n",
    " TH1 *histGlobalCorrScan=unfold.GetRhoItotal\n",
    "    (\"histGlobalCorrScan\",nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);\n",
    " /* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal(\"histCorrCoeff\",nullptr,nullptr,nullptr,kFALSE);\n",
    "\n",
    " TCanvas canvas;\n",
    " canvas.Print(\"testUnfold5.ps[\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb11415",
   "metadata": {},
   "source": [
    "========== page 1 ============\n",
    "unfolding control plots\n",
    "input, matrix, output\n",
    "tau-scan, global correlations, correlation coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b89c56ae",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.Clear();\n",
    " canvas.Divide(3,2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "951d0d3b",
   "metadata": {},
   "source": [
    "(1) all bins, compare to original MC distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "883de0ac",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(1);\n",
    " histDataReco->SetMinimum(0.0);\n",
    " histDataReco->Draw(\"E\");\n",
    " histMCReco->SetLineColor(kBlue);\n",
    " histMCReco->Draw(\"SAME HIST\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd298514",
   "metadata": {},
   "source": [
    "(2) matrix of probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea901a5",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(2);\n",
    " histProbability->Draw(\"BOX\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "994d9e08",
   "metadata": {},
   "source": [
    "(3) unfolded data, data truth, MC truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2193d337",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(3);\n",
    " gPad->SetLogy();\n",
    " histDataUnfold->Draw(\"E\");\n",
    " histDataTruth->SetLineColor(kBlue);\n",
    " histDataTruth->Draw(\"SAME HIST\");\n",
    " histMCTruth->SetLineColor(kRed);\n",
    " histMCTruth->Draw(\"SAME HIST\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ff97617",
   "metadata": {},
   "source": [
    "(4) scan of correlation vs tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0cab551",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(4);\n",
    " rhoLogTau->Draw();\n",
    " knots->Draw(\"*\");\n",
    " bestRhoLogTau->SetMarkerColor(kRed);\n",
    " bestRhoLogTau->Draw(\"*\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4758d15",
   "metadata": {},
   "source": [
    "(5) global correlation coefficients for the distributions\n",
    "used during the scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71868871",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(5);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71570efe",
   "metadata": {},
   "source": [
    "histCorrCoeff->Draw(\"BOX\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e338a82",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " histGlobalCorrScan->Draw(\"HIST\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb1a2c7a",
   "metadata": {},
   "source": [
    "(6) L-curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99c91ebc",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " canvas.cd(6);\n",
    " lCurve->Draw(\"AL\");\n",
    " bestLCurve->SetMarkerColor(kRed);\n",
    " bestLCurve->Draw(\"*\");\n",
    "\n",
    "\n",
    " canvas.Print(\"testUnfold5.ps\");\n",
    "\n",
    " canvas.Print(\"testUnfold5.ps]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d1a1563",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d78264ce",
   "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
}