{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "35603f05",
   "metadata": {},
   "source": [
    "# testUnfold3\n",
    "Simple Test program for the class TUnfoldDensity.\n",
    "\n",
    "1-dimensional unfolding with background subtraction\n",
    "\n",
    " the goal is to unfold the underlying \"true\" distribution of a variable Pt\n",
    "\n",
    " the reconstructed Pt is measured in 24 bins from 4 to 28\n",
    " the generator-level Pt is unfolded into 10 bins from 6 to 26\n",
    "  - plus underflow bin from 0 to 6\n",
    "  - plus overflow bin above 26\n",
    " there are two background sources bgr1 and bgr2\n",
    " the signal has a finite trigger efficiency at a threshold of 8 GeV\n",
    "\n",
    " one type of systematic error is studied, where the signal parameters are\n",
    " changed\n",
    "\n",
    " Finally, the unfolding is compared to a \"bin-by-bin\" correction method\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, change to use the class TUnfoldDensity\n",
    " - Version 16.1, parallel to changes in TUnfold\n",
    " - Version 16.0, parallel to changes in TUnfold\n",
    " - Version 15, simple example including background subtraction\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": 1,
   "id": "ccf6cef0",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.116209Z",
     "iopub.status.busy": "2026-05-19T20:11:23.116063Z",
     "iopub.status.idle": "2026-05-19T20:11:23.123760Z",
     "shell.execute_reply": "2026-05-19T20:11:23.123437Z"
    }
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "#include <TMath.h>\n",
    "#include <TCanvas.h>\n",
    "#include <TRandom3.h>\n",
    "#include <TF1.h>\n",
    "#include <TStyle.h>\n",
    "#include <TVector.h>\n",
    "#include <TGraph.h>\n",
    "\n",
    "#include \"TUnfoldDensity.h\"\n",
    "\n",
    "TRandom *rnd=nullptr;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60c89ef0",
   "metadata": {},
   "source": [
    " Definition of a helper function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5ff41282",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.136410Z",
     "iopub.status.busy": "2026-05-19T20:11:23.136277Z",
     "iopub.status.idle": "2026-05-19T20:11:23.143470Z",
     "shell.execute_reply": "2026-05-19T20:11:23.142535Z"
    }
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "Double_t GenerateEvent(const Double_t *parm,\n",
    "                       const Double_t *triggerParm,\n",
    "                       Double_t *intLumi,\n",
    "                       Bool_t *triggerFlag,\n",
    "                       Double_t *ptGen,Int_t *iType)\n",
    "{\n",
    "   // generate an event\n",
    "   // input:\n",
    "   //      parameters for the event generator\n",
    "   // return value:\n",
    "   //      reconstructed Pt\n",
    "   // output to pointers:\n",
    "   //      integrated luminosity\n",
    "   //      several variables only accessible on generator level\n",
    "   //\n",
    "   // the parm array defines the physical parameters\n",
    "   //  parm[0]: background source 1 fraction\n",
    "   //  parm[1]: background source 2 fraction\n",
    "   //  parm[2]: lower limit of generated Pt distribution\n",
    "   //  parm[3]: upper limit of generated Pt distribution\n",
    "   //  parm[4]: exponent for Pt distribution signal\n",
    "   //  parm[5]: exponent for Pt distribution background 1\n",
    "   //  parm[6]: exponent for Pt distribution background 2\n",
    "   //  parm[7]: resolution parameter a goes with sqrt(Pt)\n",
    "   //  parm[8]: resolution parameter b goes with Pt\n",
    "   //  triggerParm[0]: trigger threshold turn-on position\n",
    "   //  triggerParm[1]: trigger threshold turn-on width\n",
    "   //  triggerParm[2]: trigger efficiency for high pt\n",
    "   //\n",
    "   // intLumi is advanced bu 1 for each *generated* event\n",
    "   // for data, several events may be generated, until one passes the trigger\n",
    "   //\n",
    "   // some generator-level quantities are also returned:\n",
    "   //   triggerFlag: whether the event passed the trigger threshold\n",
    "   //   ptGen: the generated pt\n",
    "   //   iType: which type of process was simulated\n",
    "   //\n",
    "   // the \"triggerFlag\" also has another meaning:\n",
    "   //   if(triggerFlag==0)   only triggered events are returned\n",
    "   //\n",
    "   // Usage to generate data events\n",
    "   //      ptObs=GenerateEvent(parm,triggerParm,0,0,0)\n",
    "   //\n",
    "   // Usage to generate MC events\n",
    "   //      ptGen=GenerateEvent(parm,triggerParm,&triggerFlag,&ptGen,&iType);\n",
    "   //\n",
    "   Double_t ptObs;\n",
    "   Bool_t isTriggered=kFALSE;\n",
    "   do {\n",
    "      Int_t itype;\n",
    "      Double_t ptgen;\n",
    "      Double_t f=rnd->Rndm();\n",
    "      // decide whether this is background or signal\n",
    "      itype=0;\n",
    "      if(f<parm[0]) itype=1;\n",
    "      else if(f<parm[0]+parm[1]) itype=2;\n",
    "      // generate Pt according to distribution  pow(ptgen,a)\n",
    "      // get exponent\n",
    "      Double_t a=parm[4+itype];\n",
    "      Double_t a1=a+1.0;\n",
    "      Double_t t=rnd->Rndm();\n",
    "      if(a1 == 0.0) {\n",
    "         Double_t x0=TMath::Log(parm[2]);\n",
    "         ptgen=TMath::Exp(t*(TMath::Log(parm[3])-x0)+x0);\n",
    "      } else {\n",
    "         Double_t x0=pow(parm[2],a1);\n",
    "         ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1);\n",
    "      }\n",
    "      if(iType) *iType=itype;\n",
    "      if(ptGen) *ptGen=ptgen;\n",
    "\n",
    "      // smearing in Pt with large asymmetric tail\n",
    "      Double_t sigma=\n",
    "         TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);\n",
    "      ptObs=rnd->BreitWigner(ptgen,sigma);\n",
    "\n",
    "      // decide whether event was triggered\n",
    "      Double_t triggerProb =\n",
    "         triggerParm[2]/(1.+TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));\n",
    "      isTriggered= rnd->Rndm()<triggerProb;\n",
    "      (*intLumi) ++;\n",
    "   } while((!triggerFlag) && (!isTriggered));\n",
    "   // return trigger decision\n",
    "   if(triggerFlag) *triggerFlag=isTriggered;\n",
    "   return ptObs;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d9020f1",
   "metadata": {},
   "source": [
    "switch on histogram errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7b9644ee",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.145004Z",
     "iopub.status.busy": "2026-05-19T20:11:23.144882Z",
     "iopub.status.idle": "2026-05-19T20:11:23.459288Z",
     "shell.execute_reply": "2026-05-19T20:11:23.458775Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1::SetDefaultSumw2();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43460b2a",
   "metadata": {},
   "source": [
    "show fit result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9d892f92",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.462246Z",
     "iopub.status.busy": "2026-05-19T20:11:23.462121Z",
     "iopub.status.idle": "2026-05-19T20:11:23.665847Z",
     "shell.execute_reply": "2026-05-19T20:11:23.665471Z"
    }
   },
   "outputs": [],
   "source": [
    " gStyle->SetOptFit(1111);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b19fe9e",
   "metadata": {},
   "source": [
    "random generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e6bcfd60",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.668381Z",
     "iopub.status.busy": "2026-05-19T20:11:23.668264Z",
     "iopub.status.idle": "2026-05-19T20:11:23.871797Z",
     "shell.execute_reply": "2026-05-19T20:11:23.871145Z"
    }
   },
   "outputs": [],
   "source": [
    " rnd=new TRandom3();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a61596bb",
   "metadata": {},
   "source": [
    "data and MC luminosities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d50a5777",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:23.877560Z",
     "iopub.status.busy": "2026-05-19T20:11:23.877434Z",
     "iopub.status.idle": "2026-05-19T20:11:24.080921Z",
     "shell.execute_reply": "2026-05-19T20:11:24.080359Z"
    }
   },
   "outputs": [],
   "source": [
    " Double_t const lumiData= 30000;\n",
    " Double_t const lumiMC  =1000000;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ffe57d4",
   "metadata": {},
   "source": [
    "========================\n",
    "Step 1: define binning, book histograms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90d0e5e5",
   "metadata": {},
   "source": [
    "reconstructed pt (fine binning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b30c341f",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:24.083014Z",
     "iopub.status.busy": "2026-05-19T20:11:24.082894Z",
     "iopub.status.idle": "2026-05-19T20:11:24.286248Z",
     "shell.execute_reply": "2026-05-19T20:11:24.285737Z"
    }
   },
   "outputs": [],
   "source": [
    " Int_t const nDet=24;\n",
    " Double_t const xminDet=4.0;\n",
    " Double_t const xmaxDet=28.0;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc50f9e7",
   "metadata": {},
   "source": [
    "generated pt (coarse binning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e3fdbbc2",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:24.288073Z",
     "iopub.status.busy": "2026-05-19T20:11:24.287942Z",
     "iopub.status.idle": "2026-05-19T20:11:24.491179Z",
     "shell.execute_reply": "2026-05-19T20:11:24.490709Z"
    }
   },
   "outputs": [],
   "source": [
    " Int_t const nGen=10;\n",
    " Double_t const xminGen= 6.0;\n",
    " Double_t const xmaxGen=26.0;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e86c92c0",
   "metadata": {},
   "source": [
    "==================================================================\n",
    "book histograms\n",
    "(1) unfolding input: binning scheme is fine for detector,coarse for gen\n",
    "histUnfoldInput : reconstructed data, binning for unfolding\n",
    "histUnfoldMatrix : generated vs reconstructed distribution\n",
    "histUnfoldBgr1 : background source1, as predicted by MC\n",
    "histUnfoldBgr2 : background source2, as predicted by MC\n",
    "for systematic studies\n",
    "histUnfoldMatrixSys : generated vs reconstructed with different signal\n",
    "\n",
    "(2) histograms required for bin-by-bin method\n",
    "histDetDATAbbb : reconstructed data for bin-by-bin\n",
    "histDetMCbbb : reconstructed MC for bin-by-bin\n",
    "histDetMCBgr1bbb : reconstructed MC bgr1 for bin-by-bin\n",
    "histDetMCBgr2bbb : reconstructed MC bgr2 for bin-by-bin\n",
    "histDetMCBgrPtbbb : reconstructed MC bgr from low/high pt for bin-by-bin\n",
    "histGenMCbbb : generated MC truth\n",
    "for systematic studies\n",
    "histDetMCSysbbb : reconstructed with changed trigger\n",
    "histGenMCSysbbb : generated MC truth\n",
    "\n",
    "(3) monitoring and control\n",
    "histGenData : data truth for bias tests\n",
    "histDetMC : MC prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "595ead80",
   "metadata": {},
   "source": [
    "(1) create histograms required for unfolding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d9bf00df",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:24.493351Z",
     "iopub.status.busy": "2026-05-19T20:11:24.493231Z",
     "iopub.status.idle": "2026-05-19T20:11:24.696501Z",
     "shell.execute_reply": "2026-05-19T20:11:24.696046Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1D *histUnfoldInput=\n",
    "    new TH1D(\"unfolding input rec\",\";ptrec\",nDet,xminDet,xmaxDet);\n",
    " TH2D *histUnfoldMatrix=\n",
    "    new TH2D(\"unfolding matrix\",\";ptgen;ptrec\",\n",
    "             nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);\n",
    " TH1D *histUnfoldBgr1=\n",
    "    new TH1D(\"unfolding bgr1 rec\",\";ptrec\",nDet,xminDet,xmaxDet);\n",
    " TH1D *histUnfoldBgr2=\n",
    "    new TH1D(\"unfolding bgr2 rec\",\";ptrec\",nDet,xminDet,xmaxDet);\n",
    " TH2D *histUnfoldMatrixSys=\n",
    "    new TH2D(\"unfolding matrix sys\",\";ptgen;ptrec\",\n",
    "             nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45ccba24",
   "metadata": {},
   "source": [
    "(2) histograms required for bin-by-bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a2c09801",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:24.698395Z",
     "iopub.status.busy": "2026-05-19T20:11:24.698278Z",
     "iopub.status.idle": "2026-05-19T20:11:24.901202Z",
     "shell.execute_reply": "2026-05-19T20:11:24.900534Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1D *histBbbInput=\n",
    "    new TH1D(\"bbb input rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbSignalRec=\n",
    "    new TH1D(\"bbb signal rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbBgr1=\n",
    "    new TH1D(\"bbb bgr1 rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbBgr2=\n",
    "    new TH1D(\"bbb bgr2 rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbBgrPt=\n",
    "    new TH1D(\"bbb bgrPt rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbSignalGen=\n",
    "    new TH1D(\"bbb signal gen\",\";ptgen\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbSignalRecSys=\n",
    "    new TH1D(\"bbb signal sys rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbBgrPtSys=\n",
    "    new TH1D(\"bbb bgrPt sys rec\",\";ptrec\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histBbbSignalGenSys=\n",
    "    new TH1D(\"bbb signal sys gen\",\";ptgen\",nGen,xminGen,xmaxGen);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eaca0f1",
   "metadata": {},
   "source": [
    "(3) control histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "843bb600",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:24.902974Z",
     "iopub.status.busy": "2026-05-19T20:11:24.902855Z",
     "iopub.status.idle": "2026-05-19T20:11:25.105542Z",
     "shell.execute_reply": "2026-05-19T20:11:25.105136Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1D *histDataTruth=\n",
    "    new TH1D(\"DATA truth gen\",\";ptgen\",nGen,xminGen,xmaxGen);\n",
    " TH1D *histDetMC=\n",
    "    new TH1D(\"MC prediction rec\",\";ptrec\",nDet,xminDet,xmaxDet);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46522b92",
   "metadata": {},
   "source": [
    "==============================================================\n",
    "Step 2: generate data distributions\n",
    "\n",
    "data parameters: in real life these are unknown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1a088246",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:25.108156Z",
     "iopub.status.busy": "2026-05-19T20:11:25.108023Z",
     "iopub.status.idle": "2026-05-19T20:11:25.310640Z",
     "shell.execute_reply": "2026-05-19T20:11:25.309733Z"
    }
   },
   "outputs": [],
   "source": [
    " static Double_t parm_DATA[]={\n",
    "    0.05, // fraction of background 1 (on generator level)\n",
    "    0.05, // fraction of background 2 (on generator level)\n",
    "    3.5, // lower Pt cut (generator level)\n",
    "    100.,// upper Pt cut (generator level)\n",
    "    -3.0,// signal exponent\n",
    "    0.1, // background 1 exponent\n",
    "    -0.5, // background 2 exponent\n",
    "    0.2, // energy resolution a term\n",
    "    0.01, // energy resolution b term\n",
    " };\n",
    " static Double_t triggerParm_DATA[]={8.0, // threshold is 8 GeV\n",
    "                              4.0, // width is 4 GeV\n",
    "                              0.95 // high Pt efficiency os 95%\n",
    " };\n",
    "\n",
    " Double_t intLumi=0.0;\n",
    " while(intLumi<lumiData) {\n",
    "    Int_t iTypeGen;\n",
    "    Bool_t isTriggered;\n",
    "    Double_t ptGen;\n",
    "    Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,\n",
    "                                 &isTriggered,&ptGen,&iTypeGen);\n",
    "    if(isTriggered) {\n",
    "       // (1) histogram for unfolding\n",
    "       histUnfoldInput->Fill(ptObs);\n",
    "\n",
    "       // (2) histogram for bin-by-bin\n",
    "       histBbbInput->Fill(ptObs);\n",
    "    }\n",
    "    // (3) monitoring\n",
    "    if(iTypeGen==0) histDataTruth->Fill(ptGen);\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6a9452c",
   "metadata": {},
   "source": [
    "==============================================================\n",
    "Step 3: generate default MC distributions\n",
    "\n",
    "MC parameters\n",
    "default settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4bb6c3ad",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:25.312425Z",
     "iopub.status.busy": "2026-05-19T20:11:25.312290Z",
     "iopub.status.idle": "2026-05-19T20:11:25.514832Z",
     "shell.execute_reply": "2026-05-19T20:11:25.513749Z"
    }
   },
   "outputs": [],
   "source": [
    " static Double_t parm_MC[]={\n",
    "    0.05, // fraction of background 1 (on generator level)\n",
    "    0.05, // fraction of background 2 (on generator level)\n",
    "    3.5, // lower Pt cut (generator level)\n",
    "    100.,// upper Pt cut (generator level)\n",
    "    -4.0,// signal exponent !!! steeper than in data\n",
    "         //                     to illustrate bin-by-bin bias\n",
    "    0.1, // background 1 exponent\n",
    "    -0.5, // background 2 exponent\n",
    "    0.2, // energy resolution a term\n",
    "    0.01, // energy resolution b term\n",
    " };\n",
    " static Double_t triggerParm_MC[]={8.0, // threshold is 8 GeV\n",
    "                            4.0, // width is 4 GeV\n",
    "                            0.95 // high Pt efficiency is 95%\n",
    " };"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "900ecbc2",
   "metadata": {},
   "source": [
    "weighting factor to accommodate for the different luminosity in data and MC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "10e6bf78",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:25.516624Z",
     "iopub.status.busy": "2026-05-19T20:11:25.516448Z",
     "iopub.status.idle": "2026-05-19T20:11:25.776705Z",
     "shell.execute_reply": "2026-05-19T20:11:25.775936Z"
    }
   },
   "outputs": [],
   "source": [
    " Double_t lumiWeight = lumiData/lumiMC;\n",
    " intLumi=0.0;\n",
    " while(intLumi<lumiMC) {\n",
    "    Int_t iTypeGen;\n",
    "    Bool_t isTriggered;\n",
    "    Double_t ptGen;\n",
    "    Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,\n",
    "                                 &ptGen,&iTypeGen);\n",
    "    if(!isTriggered) ptObs=0.0;\n",
    "\n",
    "    // (1) distribution required for unfolding\n",
    "\n",
    "    if(iTypeGen==0) {\n",
    "       histUnfoldMatrix->Fill(ptGen,ptObs,lumiWeight);\n",
    "    } else if(iTypeGen==1) {\n",
    "       histUnfoldBgr1->Fill(ptObs,lumiWeight);\n",
    "    } else if(iTypeGen==2) {\n",
    "       histUnfoldBgr2->Fill(ptObs,lumiWeight);\n",
    "    }\n",
    "\n",
    "    // (2) distributions for bbb\n",
    "    if(iTypeGen==0) {\n",
    "       if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {\n",
    "          histBbbSignalRec->Fill(ptObs,lumiWeight);\n",
    "       } else {\n",
    "          histBbbBgrPt->Fill(ptObs,lumiWeight);\n",
    "       }\n",
    "       histBbbSignalGen->Fill(ptGen,lumiWeight);\n",
    "    } else if(iTypeGen==1) {\n",
    "       histBbbBgr1->Fill(ptObs,lumiWeight);\n",
    "    } else if(iTypeGen==2) {\n",
    "       histBbbBgr2->Fill(ptObs,lumiWeight);\n",
    "    }\n",
    "\n",
    "    // (3) control distribution\n",
    "    histDetMC ->Fill(ptObs,lumiWeight);\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36fd6dbd",
   "metadata": {},
   "source": [
    "==============================================================\n",
    "Step 4: generate MC distributions for systematic study\n",
    "example: study dependence on initial signal shape\n",
    "-> BGR distributions do not change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a89529ed",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:25.778351Z",
     "iopub.status.busy": "2026-05-19T20:11:25.778197Z",
     "iopub.status.idle": "2026-05-19T20:11:26.128310Z",
     "shell.execute_reply": "2026-05-19T20:11:26.127502Z"
    }
   },
   "outputs": [],
   "source": [
    " static Double_t parm_MC_SYS[]={\n",
    "    0.05, // fraction of background: unchanged\n",
    "    0.05, // fraction of background: unchanged\n",
    "    3.5, // lower Pt cut (generator level)\n",
    "    100.,// upper Pt cut (generator level)\n",
    "    -2.0, // signal exponent changed\n",
    "    0.1, // background 1 exponent\n",
    "    -0.5, // background 2 exponent\n",
    "    0.2, // energy resolution a term\n",
    "    0.01, // energy resolution b term\n",
    " };\n",
    "\n",
    " intLumi=0.0;\n",
    " while(intLumi<lumiMC) {\n",
    "    Int_t iTypeGen;\n",
    "    Bool_t isTriggered;\n",
    "    Double_t ptGen;\n",
    "    Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,\n",
    "                                 &isTriggered,&ptGen,&iTypeGen);\n",
    "    if(!isTriggered) ptObs=0.0;\n",
    "\n",
    "    // (1) for unfolding\n",
    "    if(iTypeGen==0) {\n",
    "       histUnfoldMatrixSys->Fill(ptGen,ptObs);\n",
    "    }\n",
    "\n",
    "    // (2) for bin-by-bin\n",
    "    if(iTypeGen==0) {\n",
    "       if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {\n",
    "          histBbbSignalRecSys->Fill(ptObs);\n",
    "       } else {\n",
    "          histBbbBgrPtSys->Fill(ptObs);\n",
    "       }\n",
    "       histBbbSignalGenSys->Fill(ptGen);\n",
    "    }\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a3d493",
   "metadata": {},
   "source": [
    "this method is new in version 16 of TUnfold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1d0f9104",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:26.129976Z",
     "iopub.status.busy": "2026-05-19T20:11:26.129852Z",
     "iopub.status.idle": "2026-05-19T20:11:26.333158Z",
     "shell.execute_reply": "2026-05-19T20:11:26.332523Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TUnfold version is V17.9\n"
     ]
    }
   ],
   "source": [
    " cout<<\"TUnfold version is \"<<TUnfold::GetTUnfoldVersion()<<\"\\n\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ebc3dc8",
   "metadata": {},
   "source": [
    "========================\n",
    "Step 5: unfolding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce6e4da5",
   "metadata": {},
   "source": [
    "here one has to decide about the regularisation scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db74c63",
   "metadata": {},
   "source": [
    "regularize curvature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "27b351bd",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:26.335098Z",
     "iopub.status.busy": "2026-05-19T20:11:26.334976Z",
     "iopub.status.idle": "2026-05-19T20:11:26.538773Z",
     "shell.execute_reply": "2026-05-19T20:11:26.538042Z"
    }
   },
   "outputs": [],
   "source": [
    " TUnfold::ERegMode regMode =\n",
    "    TUnfold::kRegModeCurvature;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0ed6dab",
   "metadata": {},
   "source": [
    "preserve the area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "696985a9",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:26.540801Z",
     "iopub.status.busy": "2026-05-19T20:11:26.540681Z",
     "iopub.status.idle": "2026-05-19T20:11:26.744130Z",
     "shell.execute_reply": "2026-05-19T20:11:26.743636Z"
    }
   },
   "outputs": [],
   "source": [
    " TUnfold::EConstraint constraintMode=\n",
    "    TUnfold::kEConstraintArea;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14dba1fd",
   "metadata": {},
   "source": [
    "bin content is divided by the bin width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e90762c5",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:26.746457Z",
     "iopub.status.busy": "2026-05-19T20:11:26.746312Z",
     "iopub.status.idle": "2026-05-19T20:11:26.950347Z",
     "shell.execute_reply": "2026-05-19T20:11:26.949802Z"
    }
   },
   "outputs": [],
   "source": [
    " TUnfoldDensity::EDensityMode densityFlags=\n",
    "    TUnfoldDensity::kDensityModeBinWidth;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "438826b8",
   "metadata": {},
   "source": [
    "set up matrix of migrations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "50e7cf27",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:26.951899Z",
     "iopub.status.busy": "2026-05-19T20:11:26.951786Z",
     "iopub.status.idle": "2026-05-19T20:11:27.155516Z",
     "shell.execute_reply": "2026-05-19T20:11:27.155096Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TUnfold::SetConstraint>: fConstraint=1\n",
      "Info in <TUnfold::TUnfold>: 24 input bins and 12 output bins (includes 2 underflow/overflow bins)\n",
      "Info in <TUnfoldDensity::RegularizeOneDistribution>: regularizing xaxis regMode=3 densityMode=1 axisSteering=*[UOB]\n"
     ]
    }
   ],
   "source": [
    " TUnfoldDensity unfold(histUnfoldMatrix,TUnfold::kHistMapOutputHoriz,\n",
    "                       regMode,constraintMode,densityFlags);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3bfa6bb",
   "metadata": {},
   "source": [
    "define the input vector (the measured data distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "8c861669",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:27.156728Z",
     "iopub.status.busy": "2026-05-19T20:11:27.156594Z",
     "iopub.status.idle": "2026-05-19T20:11:27.360084Z",
     "shell.execute_reply": "2026-05-19T20:11:27.359500Z"
    }
   },
   "outputs": [],
   "source": [
    " unfold.SetInput(histUnfoldInput);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a69f2cdd",
   "metadata": {},
   "source": [
    "subtract background, normalized to data luminosity\n",
    "with 10% scale error each"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "73635bf7",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:27.361329Z",
     "iopub.status.busy": "2026-05-19T20:11:27.361214Z",
     "iopub.status.idle": "2026-05-19T20:11:27.564151Z",
     "shell.execute_reply": "2026-05-19T20:11:27.563579Z"
    }
   },
   "outputs": [],
   "source": [
    " Double_t scale_bgr=1.0;\n",
    " Double_t dscale_bgr=0.1;\n",
    " unfold.SubtractBackground(histUnfoldBgr1,\"background1\",scale_bgr,dscale_bgr);\n",
    " unfold.SubtractBackground(histUnfoldBgr2,\"background2\",scale_bgr,dscale_bgr);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "138f645b",
   "metadata": {},
   "source": [
    "add systematic error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "76ffaf5c",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:27.565430Z",
     "iopub.status.busy": "2026-05-19T20:11:27.565318Z",
     "iopub.status.idle": "2026-05-19T20:11:27.768780Z",
     "shell.execute_reply": "2026-05-19T20:11:27.768226Z"
    }
   },
   "outputs": [],
   "source": [
    " unfold.AddSysError(histUnfoldMatrixSys,\"signalshape_SYS\",\n",
    "                    TUnfold::kHistMapOutputHoriz,\n",
    "                    TUnfoldSys::kSysErrModeMatrix);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "448e4a9a",
   "metadata": {},
   "source": [
    "run the unfolding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c2a1f319",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:27.769997Z",
     "iopub.status.busy": "2026-05-19T20:11:27.769886Z",
     "iopub.status.idle": "2026-05-19T20:11:27.973265Z",
     "shell.execute_reply": "2026-05-19T20:11:27.972737Z"
    }
   },
   "outputs": [],
   "source": [
    " Int_t nScan=30;\n",
    " TSpline *logTauX,*logTauY;\n",
    " TGraph *lCurve;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee266cc4",
   "metadata": {},
   "source": [
    "this method scans the parameter tau and finds the kink in the L curve\n",
    "finally, the unfolding is done for the best choice of tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "78112481",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:27.974533Z",
     "iopub.status.busy": "2026-05-19T20:11:27.974423Z",
     "iopub.status.idle": "2026-05-19T20:11:28.179058Z",
     "shell.execute_reply": "2026-05-19T20:11:28.178596Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi**2=25.3912+4.28619 / 11\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-Infinity X=1.398300 Y=5.378547\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.914692 X=1.501461 Y=5.045396\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.414692 X=1.401645 Y=5.325900\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.164692 X=1.421959 Y=5.232814\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.039692 X=1.451597 Y=5.151192\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.977192 X=1.473934 Y=5.100968\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.289692 X=1.407586 Y=5.289508\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.102192 X=1.434441 Y=5.195244\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.945942 X=1.487070 Y=5.073781\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.227192 X=1.413316 Y=5.264065\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.008442 X=1.462104 Y=5.126812\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.070942 X=1.442393 Y=5.174021\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.133442 X=1.427663 Y=5.214839\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.352192 X=1.403921 Y=5.309858\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.195942 X=1.417216 Y=5.249204\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.930317 X=1.494113 Y=5.059728\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.961567 X=1.480342 Y=5.087534\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.992817 X=1.467854 Y=5.114066\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.024067 X=1.456686 Y=5.139191\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.258442 X=1.410143 Y=5.277471\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.055317 X=1.446834 Y=5.162805\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.086567 X=1.438265 Y=5.184835\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.117817 X=1.430912 Y=5.205245\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.320942 X=1.405542 Y=5.300271\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.149067 X=1.424684 Y=5.224027\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.180317 X=1.419475 Y=5.241204\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.383442 X=1.402644 Y=5.318368\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.922504 X=1.497750 Y=5.052595\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-1.938129 X=1.490553 Y=5.066790\n",
      "Info in <TUnfoldDensity::ScanLcurve>: logtau=-2.211567 X=1.415168 Y=5.256821\n",
      "Info in <TUnfoldDensity::ScanLcurve>: Result logtau=-2.336567 X=1.404684 Y=5.305205\n"
     ]
    }
   ],
   "source": [
    " Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);\n",
    "\n",
    " cout<<\"chi**2=\"<<unfold.GetChi2A()<<\"+\"<<unfold.GetChi2L()\n",
    "     <<\" / \"<<unfold.GetNdf()<<\"\\n\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "374660da",
   "metadata": {},
   "source": [
    "save graphs with one point to visualize best choice of tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "45ebfeec",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:28.180212Z",
     "iopub.status.busy": "2026-05-19T20:11:28.180103Z",
     "iopub.status.idle": "2026-05-19T20:11:28.383564Z",
     "shell.execute_reply": "2026-05-19T20:11:28.383014Z"
    }
   },
   "outputs": [],
   "source": [
    " Double_t t[1],x[1],y[1];\n",
    " logTauX->GetKnot(iBest,t[0],x[0]);\n",
    " logTauY->GetKnot(iBest,t[0],y[0]);\n",
    " TGraph *bestLcurve=new TGraph(1,x,y);\n",
    " TGraph *bestLogTauLogChi2=new TGraph(1,t,x);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ff533a3",
   "metadata": {},
   "source": [
    "===========================\n",
    "Step 6: retrieve unfolding results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9505a3d0",
   "metadata": {},
   "source": [
    "get unfolding output\n",
    "includes the statistical and background errors\n",
    "but not the other systematic uncertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "8477af89",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:28.384897Z",
     "iopub.status.busy": "2026-05-19T20:11:28.384783Z",
     "iopub.status.idle": "2026-05-19T20:11:28.587628Z",
     "shell.execute_reply": "2026-05-19T20:11:28.587062Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1 *histUnfoldOutput=unfold.GetOutput(\"PT(unfold,stat+bgrerr)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "294c618c",
   "metadata": {},
   "source": [
    "retrieve error matrix of statistical errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "93db06e2",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:28.588900Z",
     "iopub.status.busy": "2026-05-19T20:11:28.588788Z",
     "iopub.status.idle": "2026-05-19T20:11:28.791683Z",
     "shell.execute_reply": "2026-05-19T20:11:28.791106Z"
    }
   },
   "outputs": [],
   "source": [
    " TH2 *histEmatStat=unfold.GetEmatrixInput(\"unfolding stat error matrix\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efd6e475",
   "metadata": {},
   "source": [
    "retrieve full error matrix\n",
    "This includes all systematic errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "799905e5",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:28.792924Z",
     "iopub.status.busy": "2026-05-19T20:11:28.792810Z",
     "iopub.status.idle": "2026-05-19T20:11:28.995628Z",
     "shell.execute_reply": "2026-05-19T20:11:28.995103Z"
    }
   },
   "outputs": [],
   "source": [
    " TH2 *histEmatTotal=unfold.GetEmatrixTotal(\"unfolding total error matrix\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1821b7f7",
   "metadata": {},
   "source": [
    "create two copies of the unfolded data, one with statistical errors\n",
    "the other with total errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "e915cae6",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:28.996907Z",
     "iopub.status.busy": "2026-05-19T20:11:28.996793Z",
     "iopub.status.idle": "2026-05-19T20:11:29.210249Z",
     "shell.execute_reply": "2026-05-19T20:11:29.209669Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1 *histUnfoldStat=new TH1D(\"PT(unfold,staterr)\",\";Pt(gen)\",\n",
    "                               nGen,xminGen,xmaxGen);\n",
    " TH1 *histUnfoldTotal=new TH1D(\"PT(unfold,totalerr)\",\";Pt(gen)\",\n",
    "                                nGen,xminGen,xmaxGen);\n",
    " for(Int_t i=0;i<nGen+2;i++) {\n",
    "    Double_t c=histUnfoldOutput->GetBinContent(i);\n",
    "\n",
    "    // histogram with unfolded data and stat errors\n",
    "    histUnfoldStat->SetBinContent(i,c);\n",
    "    histUnfoldStat->SetBinError\n",
    "       (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i)));\n",
    "\n",
    "    // histogram with unfolded data and total errors\n",
    "    histUnfoldTotal->SetBinContent(i,c);\n",
    "    histUnfoldTotal->SetBinError\n",
    "       (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i)));\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93e3dbaf",
   "metadata": {},
   "source": [
    "create histogram with correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9c9df8c3",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:29.214174Z",
     "iopub.status.busy": "2026-05-19T20:11:29.214052Z",
     "iopub.status.idle": "2026-05-19T20:11:29.418891Z",
     "shell.execute_reply": "2026-05-19T20:11:29.417770Z"
    }
   },
   "outputs": [],
   "source": [
    " TH2D *histCorr=new TH2D(\"Corr(total)\",\";Pt(gen);Pt(gen)\",\n",
    "                         nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);\n",
    " for(Int_t i=0;i<nGen+2;i++) {\n",
    "    Double_t ei,ej;\n",
    "    ei=TMath::Sqrt(histEmatTotal->GetBinContent(i,i));\n",
    "    if(ei<=0.0) continue;\n",
    "    for(Int_t j=0;j<nGen+2;j++) {\n",
    "       ej=TMath::Sqrt(histEmatTotal->GetBinContent(j,j));\n",
    "       if(ej<=0.0) continue;\n",
    "       histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej);\n",
    "    }\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d753a09b",
   "metadata": {},
   "source": [
    "retrieve bgr source 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "d5f0d591",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:29.420195Z",
     "iopub.status.busy": "2026-05-19T20:11:29.420081Z",
     "iopub.status.idle": "2026-05-19T20:11:29.639301Z",
     "shell.execute_reply": "2026-05-19T20:11:29.638659Z"
    }
   },
   "outputs": [],
   "source": [
    " TH1 *histDetNormBgr1=unfold.GetBackground(\"bgr1 normalized\",\n",
    "                                            \"background1\");\n",
    " TH1 *histDetNormBgrTotal=unfold.GetBackground(\"bgr total normalized\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93fa6c09",
   "metadata": {},
   "source": [
    "========================\n",
    "Step 7: plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "95379bb1",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:29.650233Z",
     "iopub.status.busy": "2026-05-19T20:11:29.650093Z",
     "iopub.status.idle": "2026-05-19T20:11:29.887021Z",
     "shell.execute_reply": "2026-05-19T20:11:29.882387Z"
    }
   },
   "outputs": [],
   "source": [
    " TCanvas output;\n",
    " output.Divide(3,2);\n",
    " output.cd(1);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "433198b6",
   "metadata": {},
   "source": [
    "data, MC prediction, background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "f2791517",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:29.900937Z",
     "iopub.status.busy": "2026-05-19T20:11:29.900800Z",
     "iopub.status.idle": "2026-05-19T20:11:30.104599Z",
     "shell.execute_reply": "2026-05-19T20:11:30.103965Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldInput->SetMinimum(0.0);\n",
    " histUnfoldInput->Draw(\"E\");\n",
    " histDetMC->SetMinimum(0.0);\n",
    " histDetMC->SetLineColor(kBlue);\n",
    " histDetNormBgrTotal->SetLineColor(kRed);\n",
    " histDetNormBgr1->SetLineColor(kCyan);\n",
    " histDetMC->Draw(\"SAME HIST\");\n",
    " histDetNormBgr1->Draw(\"SAME HIST\");\n",
    " histDetNormBgrTotal->Draw(\"SAME HIST\");\n",
    "\n",
    " output.cd(2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3633e1",
   "metadata": {},
   "source": [
    "unfolded data, data truth, MC truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "6223a9a2",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:30.105942Z",
     "iopub.status.busy": "2026-05-19T20:11:30.105825Z",
     "iopub.status.idle": "2026-05-19T20:11:30.309418Z",
     "shell.execute_reply": "2026-05-19T20:11:30.308855Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldTotal->SetMinimum(0.0);\n",
    " histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53c9862f",
   "metadata": {},
   "source": [
    "outer error: total error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "3ba809bf",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:30.310788Z",
     "iopub.status.busy": "2026-05-19T20:11:30.310661Z",
     "iopub.status.idle": "2026-05-19T20:11:30.513566Z",
     "shell.execute_reply": "2026-05-19T20:11:30.513020Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldTotal->Draw(\"E\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ef3419d",
   "metadata": {},
   "source": [
    "middle error: stat+bgr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "4d34b4d6",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:30.514932Z",
     "iopub.status.busy": "2026-05-19T20:11:30.514820Z",
     "iopub.status.idle": "2026-05-19T20:11:30.723437Z",
     "shell.execute_reply": "2026-05-19T20:11:30.722725Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldOutput->Draw(\"SAME E1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c449ecd1",
   "metadata": {},
   "source": [
    "inner error: stat only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "6764d006",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:30.725385Z",
     "iopub.status.busy": "2026-05-19T20:11:30.725259Z",
     "iopub.status.idle": "2026-05-19T20:11:30.934790Z",
     "shell.execute_reply": "2026-05-19T20:11:30.934080Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldStat->Draw(\"SAME E1\");\n",
    " histDataTruth->Draw(\"SAME HIST\");\n",
    " histBbbSignalGen->SetLineColor(kBlue);\n",
    " histBbbSignalGen->Draw(\"SAME HIST\");\n",
    "\n",
    " output.cd(3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14708b2f",
   "metadata": {},
   "source": [
    "unfolding matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "75630570",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:30.936714Z",
     "iopub.status.busy": "2026-05-19T20:11:30.936555Z",
     "iopub.status.idle": "2026-05-19T20:11:31.140287Z",
     "shell.execute_reply": "2026-05-19T20:11:31.139074Z"
    }
   },
   "outputs": [],
   "source": [
    " histUnfoldMatrix->SetLineColor(kBlue);\n",
    " histUnfoldMatrix->Draw(\"BOX\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "323a7982",
   "metadata": {},
   "source": [
    "show tau as a function of chi**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "8b9244c4",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:31.142040Z",
     "iopub.status.busy": "2026-05-19T20:11:31.141899Z",
     "iopub.status.idle": "2026-05-19T20:11:31.345754Z",
     "shell.execute_reply": "2026-05-19T20:11:31.344629Z"
    }
   },
   "outputs": [],
   "source": [
    " output.cd(4);\n",
    " logTauX->Draw();\n",
    " bestLogTauLogChi2->SetMarkerColor(kRed);\n",
    " bestLogTauLogChi2->Draw(\"*\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1dc173c",
   "metadata": {},
   "source": [
    "show the L curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "84eb261f",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:31.347286Z",
     "iopub.status.busy": "2026-05-19T20:11:31.347153Z",
     "iopub.status.idle": "2026-05-19T20:11:31.550575Z",
     "shell.execute_reply": "2026-05-19T20:11:31.549551Z"
    }
   },
   "outputs": [],
   "source": [
    " output.cd(5);\n",
    " lCurve->Draw(\"AL\");\n",
    " bestLcurve->SetMarkerColor(kRed);\n",
    " bestLcurve->Draw(\"*\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d2cbd0",
   "metadata": {},
   "source": [
    "show correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "0d97996a",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:31.552030Z",
     "iopub.status.busy": "2026-05-19T20:11:31.551912Z",
     "iopub.status.idle": "2026-05-19T20:11:31.755837Z",
     "shell.execute_reply": "2026-05-19T20:11:31.754988Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "<div id=\"root_plot_1779221491743\" style=\"width: 700px; height: 500px; position: relative\">\n",
       "</div>\n",
       "\n",
       "</div>\n",
       "<script>\n",
       "   function process_root_plot_1779221491743() {\n",
       "      function execCode(Core) {\n",
       "         Core.settings.HandleKeys = false;\n",
       "         \n",
       "Core.unzipJSON(88846,'WkwIy0wADlsBeAHtvWmTHDeSoP1XaGmzZr37BsPguJFh84GipFG/q4Omo0Vu29hYkpUkc1Ws4mQlu6ke6/++9jgQkZFZSTUlHiJ7KrpZSiAOOACHw+Hnfy3+Y/fz8/XF6tl6sVx8f3d18ZfV1Y/rh99drJ5fPb3cLbrF4x8uNv/5Yv3HTxdL0y0ef7LZXdVf3zz8v+tHO+oXPPbN893m8qIV/vfm4myxdN3i8fSl5X+dautVDTjnrfO+Wzz+cnOxvnt5frldLKUVv9v9fL7eF3/cnO2e1uLnm/Pz9jDAUhwfNkZfXz/efbXaPtlcLJamp+bbzZOnR1WfXO52l88OH/v+8vlhxf3HG4Cw3eLxg/3PO/UnH75/tVvtaKUUnjko3aklXv58u3q2PoabuqOOT88ddmh6dKweP8nr85Hhk59cbs/W2+82f2ujN6v86vJsXef1vtT/Pmj/vW/r2D5o/72/u7zz8Ore5uX6/Kf2xu7yVHGxjCXS9YMXfNIh279yWLFY3q4VPxy8xWz+sH9nXhyb+dPBC/Ujf9q/clgxNTN2ZXd5v/VGu3K9yEwaI97FlCWaIiWAnvv3H/xUx6m9TxEwx+JieZsPWJEcvTcxpRKj4sn55V+//vRuG/R54f4Pz/UGTz2Y/f5xqv1i+nXn4dXBh+48vDr41p2HV/vX7jy82r/5w8tnuhoY4Z/3P18+W72sHfrh5+nn90/Xu9Vi6bRjTzft152r5+tHu29Xu81l7cXXL549XG/r7+83j356uf/5c/355eWTVvnl5ZN93d/q3Xurs3urzQWrp1s8vru9vLp6utq0D07Fe5eNEs3RGuSq5T1Gf3V5tnm8WZ8tlo9X51frbvH437abs5eHxZ/3xTsPr+5eXm5nz392ttmtHrLed9sXfODzzcv12UG/x0/f226ebXabv6yvrpG8LzdXUNSR2rbiartdLP/8793i8vmOH3/vFo8/e7l+dLVYXrw4P+8Wj7+u5PmRQC2/3+yAo5W+fvHs3up8vduNNJIB+3r9cne99tM/fnfvyzsPFsvFv4w/u8XjTy9fPDxff/Li8eNxwr5d71abC0ar9fX+1eZv6x+uxvsPDot699v16nyxtDSut2tZfG9zSKDvj5uLs8u/fn/5nHU0Lz+Ylxu92j/wxRrq3NDgr+NKv/t00Zbz3dVud22Y7+x2dSOjf/c/We/+ul5fNEp9UNKh/Hx7+ez7y+eLpfTgzv2z1Q5CqIUHY4Gd7E4tyN+7xU9fXf5l/c3z1X++mDDip2/XjMth5eMvNk+efkkX2qak2LnaPXo6Du5P3z29/Otnf1lf7L7brXYvriYk/OnOi90laDA9+dX64sUnq20tgyR3HoFm0xuPv12vzr65OP95fOPxj5vd08sXuzlGjlj6xeqq4dhYM3/qz0d79VvjCKCYr+QIflw/1FW/uXjyKrYAzLh7vrq6akuC5yofMq94DqouzNKG0LV/gyxNZzozWK3ll1savWsGr/dsCEOY3jFDbE9Sn9qz/M7LXDor0mU/FH5n31mJg9Cg79q/QWQpxXaSbSfJDmKXKXb6/0HcUozt2r9B/FJs6tq/QcJSguvav0HiUlLu2r9B0tIa37V/g+SltaVr/wYpS+tt1/4N1ujDUnjfDFYOi3YpxXVSpJOYB+uWkmMnmaIbrK/txtBJMIMNS8mOO52EMFigakB6M9i0lBA68bYTZwabl+KCdjDHwZalpNBZE4FkcGYpznbadHSDk6W4BASdRDM4uxRnOgm548PO1XZyru8yVLa2E9zgwlIkdfq+D4OLSzEeCBRkl5ZiSqcvODu4vBTrdVytKYMr+1EPYfBmKcl02t0ig5clU6wwuzB4u7RWOh0figxVrkA6M3i/1CFm1qwffFjqmAKGlMHHJZ8pvitx8Kk2w0D5NPi81BZFwRx8WYpPnXZC7BCMAqEoNgRZio+d4uwQrI5DK7hlyeNvv0x+/B2WQcbfcWnT+Dst3fgzg9jtM2VpOg8ChyGapelirr9ZOkXqb7ukU7qshsgKElfqHRYR2KGvBwo51UKkDamNRF1JrrWiresilSHSvBZsGhLta8H4IQEABclmSECghRCHNK5h5jcBgd4xeUhAQCGHISkAIXRRhjSuZJeGNDYvbkhlWSHszJDNsgLIb1kG7Qa/7TLrMPCbJVypixkyWKnDQKEuX9ozQ2aZxLHA0q1QmSGzcHUYKLBsdbQ7M5RGuFymIJUiiadQiZcUyFdxekdipOBrwVsKlYqBVGYolYwV/XLSp5JQn/W3T/wu+tu6zgxizFJEOh+UKIgRxVBvO28HMbYShcL0yiDGLXPqRLGG2177ISF2jmJYWps7l7rgBzFxKcF0iQVaBjGQjKjPSuRTLE/bOddZS5Eh6SKjKoMobQ2sITuIQMakc9JZSlZHFWIUyyDQVn7zz4ZBIK7MJgPPdyVUiCs+C9TVSFdKJ8YMImkJgaM/kgeBuoJ3sXQCRS9LyYWGZRBrluKhWb6TCH2HjtVxyxTtUl8zqctlEOuW1mk7ZhBIaw5sCx3LXKCt4jolZZaHoy51K7nCDHWdbQHWyLIIS7LIYJkSSGL9N1gDbXJd+zdY5kToj/4brAmV3kLWB8ukZKckwAzWpCVrM0mX+HJeWlcqWJIGa8qStyDWYbBMCVjAy4MVWSbp6suDFat0ub09WCaF9nh5sOJ5s748WAm8WV8eLPMBRZfE24OVVKEdG81T57TVovCOzbLdQY4q0IPVKQHa2i5TUvsq2Q1Wp0Q7a50frPVTb7VlJqX1FqAtVHw/qjaN/a0v56nDvOzGNQxtM35w4zIGmdi3WMltF6/lupj1vrGDYz2P97U84wm0XFc1zwcZHAu7Pa7FPHEQWoS2VXYjyODNnjEygzeVJ+B1SnZ6k9KeT6Hk2U7YUQZvGI1Kveif18WstF/741nMlNmNgxs8qGPs1D8vppa5b+zgYZYAkR1Uy5a2Wve8OC0Zr/CLwsGj9EZ0k+NlLY2k3gxelNTXfgmknleoZ6PhBTN4yz5Tu2TblFSecfCwQ61LbGmeRdu6pF0GP1qX6n0oidUu1/twH61LDIHNY5fq42Xskz7tDEUa0IedaLHO2uDduOnpq053vTrfg3e68U4v6s47vcdw6KQNAZRsOMJXQpFpqmt5P/VAFMp+8ut9X/tT+zgEULK2oyCHEgG5ztoQShpL+i3te+3gEIp2nZKxQzRMge6KQzSVS+edaOizTtoQDT3WSRuiqdw6HANYWFeYghghX3WFtXLWLrdpGyJ7CvfrtA1RSZhiYG0SLKwYqF2KbC26wurnxOnn6CIAsrMoEraH63iM3ZL58SE2rh2qx5t1dbTOsauAjXU4rCICszZEUHLsH+9BsVr/wJIIyaorrJXbxt+mPyrNqv3ToVYuXVdYK9ftrSJItGXsHm3Dpdfeaamuz/ak2x+ioquDAkcS4c5ZY/q7Dkf9rajBzzQSwgr2kBohbMM+JLCugRhkSCBdK/JymPqrN/fTT/dSY28UhCE1BoeW9OHau3ozjzRQqUOGBrZxDTJkUymBTsqQGxFUIpHNNOfcaeBU2pMbNjZ4hzxio1KaDC423NM2psGGx6wUuTbRkJCHeZCttJEVnlRK2D65J4TcaeSodjezj44E2dghs49WtNeH81TSNiZaBJOrCEjzFNruUEdnyGyhI13nsxBJsK0S4wyNbEVebpxwA6lhYxuFITdsHFuqa5USIHFmbA3BbDfaVYECGVvfWBaZI2MbIpZFdnWQKt4MeU8ggdjpfgFM9dlGtdtGnd14PqgP72m2NqRHlPFVPz+jZK8Um+/od71S7Lr4huz1pNTwxI8Em376PThA5/fg1M9M4NAxP4FTHx7B4UNhhEYLbYTqFj1kjottXrQjYUIj/W4Aokadcz01tgEdcqgQtQEfcgAivYsoJChE+yIQNfqY9QzZ9oMc64ZWp23I9RxZJ23IepBshDPXg2SFdsixbu8VFYYcZ5ss7etxcv/VOj5jr/VEOb1aj5TjvdREISNAqe1+I0ipIvY4KKkidqXYOVXEHjuaKjM0jULiZDOesHOqiD0OBAfNehOMyhWx6/aXc5228dFcZ61uGjnrnNVdImflgRo4nDgVt+rkcuZsNFPHKNelP7ZRl/7Uhi791kRpe4++1ij1OAKNULf2G3Og21RuRLoVKlEcZ7PR6Gk+G40ePzNtQIrQjUbrh0oj0a2wp9CMW5mRaDC4QKTbXmzsUBqZrq2URqdbT8pIp9tUl0an62CWxjOMpcoxtJkujVC3m026pjgFUBw3GPj6bygQS/C4/huK1Rmse2ixfhkcx2pxvKpo7jkap6HYuAxGT+CeQ39a+sJBXTwfzUufOw6ePg7FlqWPXTCd+DwUZ5Ye8tlJkKE4WXrfKTB+KM4uveuUPKWhOLf0tuP4HZEk+KU3HRK3aIfiwtKVLhRkYUNxcelyF6uIsLi0dKnjiJxoIS9d6DgAJ1ooiE1i7CTFoXizdK6LqZNUhuJl6Synec7Gxdul0yOf5DAU75a2dJyFcxqK90ubVUhQEGsExFYpIC0cCpK02KXUSaGFhIQglU4KLXBY7TKYLEPxZWltlytfV1SO1unRLw0F4li6HJBJDQXSmLscOytuKAHpQZdzZyUMJfilxK4gvcpDCXogRBBmZSjQxdAV20oseRVPWL4Jf4fIpZb02FxiLUX4qq6kVgJlusI5H1ENMgNkEVbKUKKjJAgutKjAKNZqMdS7rt2NtVgFFyWmWqSTeSgx0y0x9JIiAtkO0YQWE6LPDlGplTSUJFpEHKZFRD4dIqFa1BESsa2IlLhDeoEIvEAZcycqa6CIdKHKC7WYapF2ucthvhOkHIw2lJEiYIShZKNFJAValFoEDO4iu1ahRy26WgQMPxQEctyl+xQRyHZi6T7FWIuAQTHVImC4oUAgeRgwKBYtOsBwQykVKmRPYofCISojfWzFOlbIm0SGApVEOg0YFOsMOsCgWGfQAQYyOEUocYBBUTFKpataVJQSD36XoXCk950gxTLI0hDaUQbhM2WEURCMsYxUvROPOELvu6UggY6dNYnnkWAikR7LoZbzWEYs04kvSMR5Pmk5AI6W2VUgO501gfsqTpUAPJTFLEvbirUoWlQ1BnftEmF6lZWIEadFYNG7XouA4rgbloh0gESLkWIEECSOCPCQJ47FrEXA0LtliTSPUUHEaI2+q9JBikJDkLPCmFoLkKg/CkOmrK6KW6UwZPC6tpXpE9wuMkioE2CrEK+VAVTZ3armKMCCyAjIoV8KTEV/pZE0zwGsoLbpJNO+g2NpZdrnEIZqSbUSlJEwVpliBh6HLLaVgcdxOKxyUeSTxqlGQOWkVZKbILyCeDYzGi5DllHcQMHFQN5pnzLwQOBpv9J7MZB42qdM+xB52mcHoH3v2CyqHJPx8J6tpJaBxwc2GoGeJ9r3kW1IMmNO+z6xSdWytp/Zwqq2ivHwhQ2ulhmPYNj+VP0WgScIm6OKRSPtB8tGWsu033hi1XExH8EvA++zodM+fHHqJHNwoP0Ql+yRWtb20jICD/e1vbwEd7RM/0NZouJC7qui8miQINYy7UVZIrPmvqe9qArDWqa96JaMDRo6T3vRL1Etapn2YljSlpZpL0YVdGuZ/sW0VKk9fIS2l5co/yirOD6WZUH7R5n2kkGrVMuML9uCAWDoHQ2yMXDA0QpaVJ6ZIVKKSIVfso3UCtpke2DR5KJqQjFsEEhZqbDaKoSFZqnQZlEL0CwV2ix6O9aoUWWbGDYKi2qRCmaWrQLU1QqaZbNQjSNPgFtsFyAjmk/k/IYNA2zkCaFZtgzQb1+BTB6tJU8w/mwbIJhWABgbh29woC0wbB2gmD4BYGwecGdawQCNMjitAFIVwrXxECBVeUgbD6X0uokwpmgsAF23kTamSsvZSBTzeAJI2UogS1oBpGwmlc6gaqQChqDNLRWqBFJKkDuBzKoaiKWtFWhdRk2QVqiGBr3D+A1UKewpqojIHfQUbZBiM2V0K6p6oGu5g6CqQoiByAlFqahCSBXoqYPAohGSAr6lDgKLTkgK3UodBBa1kKDHoaxaIsQ4ADSWq6RKywCsagi6BItHmZML8ETV2KAaQmvJ8xBgdEOohLUMfCiHFJkzSkZBO4SShf5AcFEQWcXt3GXgQ5YCxUedS/tI9xqqQ3BVRQSFL6aD4KIlUp0zZdpHugfFL9Ilxgt5ii4E6SC4aIksFJ6yto+GEZy2HQRXkO9B4SkzPghVdJm4DoIrbDCotorrqiLNLq2uGt9BcEU3GBDedxBcYYPRRQTTTzksrQP/ORJQjkt0mFqmfTYYXVKxg+AKG4ynvYhiVsRV2wXmFIIr3iytLrDUQXAFlTwUnDLto5PX9ZY7CLBwjtDlljsIsHCS8Ho66ALtc5bQ1Vc6CLBwmlAGRY+YIpwnVM5gusB4cKLQsu0gyMKZQsuu87TPqULLoYMgS+CMCmsROwiyBKSkyoB1EGQJSG/hRErngSegs4VxFFWtSgA+yq6DQAsqey37DgItAfi4HzsItATgo5xQZosE4KNcOgc8Efg4rZgOei0R+CjbDnItEfgo+84BTwQ+yqGDWEsEPsqpc4xHrOMHDqGklAh8prOwxsATgY+ydFBuicBH2XVW4QE+yqjNUbUCH+XYQbYFkYyWOQ1RBj7uw4lTZn5LB85AxCXV+cXWRYAnhXafNUg5tjIcCWXOi7wP10UZ/KOMchdDok9ePMaW7F/UIG/x+fnlaufsolucq/VYCN3iL4vln4v1XbEc/WJXLLQld8WWrjjTFVS8znbFua443xUXuuJiVxw0KnfFla540xUvXfG2K951xfuueA6P2IhA23JXfOlKMF0J0pVguxJch7lpQQceYlcCNDF3JZSuRNOVKF2JtivRdSViaBK6EmNXIrQ0dyWWriTTlSRdSbYryXUl+a4kjqKxKyl1JUGHS1ey6UqWrmTbley6kn1XMofU2BVoqdLv0pViuoIOutiuFNcVDFxK6ErhCJu6UqD7qldHDgMlNioJMWggVDJq0MVAgg1sg4H4Ggx2DGTXQGsNBNbAVxtIq4GeGj0VQ0kN5NPAJRsIp4FPMJBIAydsII4GimhgDwy00EAADVTPwMwa6J2ByBl4VwN5M2q4BCEzcAQGEmagWwb200CxDChn4DYNyGagSgY0M9V4hzfAKoNgxkBzTOANlcrBVBr2eIP4xcBTGnZ4Azdp2NoN0hbDnm6QtBhOJIYd3SBcMWzlBlbSsIcbDgiGzdsgQjHwkIat28CpG/ZsAw9uMCAwcJAGltqwVRsYSMMebWAdDbyngWE0cIFGt3ZlGtiWDcySYUM28I2GrdgU3oBtNOzEBobRMOdILWAa+IORA3NerZqYcwQR7Pf84cTPnCujiJRBlEOs9hDMufKHwpxXwwDmXO2lEA2wq/JHDSlogznn3M/2yB/aYM6VIawWWcy5soPCnHNYZ+/iD28w5xzC2Y8QZNAGc67MHwdqNhb+8AZzrgZbyudxPGZb4A9vMOcqDBTmXO3j1OJFmHM1cOMYC6HmD28w52php6fQahnGnCv7phZyyrfpoVAZtmpCwpwru6anM+XT9NwlKf/73//+9+5dWXOGX7LmrP4Iv+A/0iw2MdPeXvDgaNlca25NHzhyhKA4c3Y49nPYe4AIfhfXXUC+Wm1/Wm9nLiW1YvbJVjG5SXy/frm7c/EEs2sMUCnWm6Y3jIHeP988uVggFanl2fe5/fkl1uxRLYxXLzfXbcPv7HZ3qMd8+2zzl83V5vLiarEMQovcmX3wy9XD9ejmQntari14WtDyN48fX63V/wQ62yonsJ3CvXn005friyd4z5jeYOSsczC+qn3B3Pz4td35aG4+PTI2j2Hyg4+mhzqdv6GH/+ej6eE0Qb9yDj9ZbWfeQ5+stiNSqAk2/lgs0/N739U18el29dfqdlHL3zzf7V08aqF5edRCc/T45vnu02pvXz3IMIJnEeky+ub5rpEEOvHN893n6nHVHv18s3/ymhE9z2jl2WaH59hY/v7y8lxt6KmoXip3Ly92ly+2V8194c6uQXRENO/sdqxjpVO/QA7sr6QHLBdGoLoX0U9KOEX0Wvrs4uyz7fay+WyxuLWoj9PU5y8uHjXKwE2KM0JGsc0id/F6aQ/Tf4rtYZY+xdmUf7l+sr44m7vWAF2tnRFZPrSvHNsene34xEgZpgcZv4aN3eLxF7hDrK+OyHir/e756hHuANr25AE368Pk/tbqgHF67hCa6dGxenz0qGl97rjfs8q9X9EXmytwcg4PVXyvgRMN3R6fGxuuozM+Ota2B4+g4amvNhebZy+e/Z/19nLv7cGNA6dEpfLV8eXedv14vf23L/dP1/rZwNWKeTeBdF6772et/XT9+IvFMhgmbar5cbFMhzX3F6gCZo88aBX3VnPku7c6wC0an6r2LWvVdWfMe6uzg64zdvdWZ9edO++tzk74d95bnYHs9/fD02oeHNSwPTZ3JRrcPPqpOSvdWz2vXpT3G9GYKh4s0FAuHn/3aLteX3y+eqTUB/CgbLPhp8g6mKEtVfP5GN+arR8eobhfPdSMCFTRqj6zfQYhXcTeqy8NlUozSm6wKInhfFpB+3F0eeLJLyiI9mO33Tz/dP1o82x1fjX5FSlVbgyPnbiGWe/0gaPuad28fyCJVs46OJWn3aryKfAg1el3emsigfpO7U6ZnqY/4/fpTiU+l3+8uFhvv6V7PMlS089eLZZ/xpfo1i3+iL0lnPtv+Vt5LMdb/pZEStbwgG+P7EtaNa+35vDx6V427V42t7y55W95fVLbfkd//h3isl6drbds2Oo2pcM2lT7f7D4fkSY0pFFnLGZxuqO4pKP2aHWuLzP7///l5oLKkRe4u3o+L36/eTZxlCnnItkrxfjjs9WTNR+aCPzd1cXZ+frHp5urn9bbb1cXT5qbcq3/5PJlq6uzV2sVkpmn5p82l+ebi7G2eTDWR+9uto/Oj6l9u4WLKUDPNsD7cN2fvXx+f875jJUP5pUPTj05Vh48yYNfrV5+unmibvQg4Tfb3dPLu6tn6+2qUZ/rJ7Z7q7O35oLHjL3SBe/e6uxVRzYOU4zQjGJRHIlPuzsbv6Od8fpB7Dq5PkGrb3zxCSDAlVL02Yao0QRksbxtg+uRO5WQI+oIvFftYumkN3pZ6wJigYW67lubUl+85cmcChN26Ml/2/TOFOckxSg2e+j/zLP/dupt5HZEcoWobH97sWw3nXXBZadHyUO3f2tyH60rJSfnkA9qXIAxbIA1ri+Gu9kVf3C3+s3nEsS6ZJLLFvmORgMYfe9/WixjXwKDlL33DhXML4QLsJJ7kxkih5yJfh7GDrBOeqsPqKROb0+Q2tiHyMsI7ve3mI+jWyN4LZZAb7wNwZYUfRIfdWdrwQE09MCJaR5v637sUoheSvYpx3I99IB1zvQxWRusz5lV3t4mEMEpXBlvL5a3RXpJEr2xJnnkW4vH+2ACvTEVAUcMZHD3sQX6cHCbhse4BaZ31tUrhJgD4zXGMaC7FU0bnuIUXsMg9E7aWykbF2muRiswvddZMcZ4n6oM4yDuwUlQDyIhnIB2ioxwquEpUsKptlvkBIbreuQEqzhYQydYY32fXfISjY2JFfwrIykwdDeRFKZANb86ksLRaf71359iL/zHUfSFWn7d+Atvz1+/sTTvxV0frHslr6BH4ldxCyD/2+UWCEvDSpuHpdE1RlyaEwtsftSgH7XcjpbXWaw35K8+YyBq1KNfHLQv5NNfGrKGbS8uHl+en20untzaXDx/sbu1XT/iLWXZG/s+G9wqZHhDZuzXS8W/frQ+51CoB95XCLWbRPtVIaRad1+qNHXWwee71uOPWhSOiPHhBim+1R1VY+vUX4TWUeS9Xx8g/Mvjzzfbqxbr5MvV+It4W3Y8zD5bf7q5en6+msWq4XQ1nc8YQxX+72PHfHV59uXqYSv/gnD+9Wbq5+OZosUPZpKmE/8kPWhqjl/QV+wniXVb54jv3J9iHr3/KXqFduH1puhvH/QUMchKxj7yKTpSiMy0JUaFG59d7LYb5Epik4Mn/v7qxbO/LpY5CxRAS/aw+HKxzI6jhws2+ig+K5/Iey8JvWaiK8H3ko3lX7d4/NXqJXJhZeD1MNIExZVcfH25fTZKIMHopuiosaYef/fiGQAcmrtE7GCquUus1i7i1QgUywCLgUfA/8J3IRMqBMOj2DmMnAjUoFb/mPRipKVRETSSgmr2rcVzQk3rCIKTDa4oMWC9Ixrx6iB8IHoLogleV5P+EuN2xOPdW/1ljTTzlTEG375EAx7F+t7rwSaaVGI9fstiKd74PjlTctDDLCSGw7ppT6dMXDkwRfkZEdd7600KXrJG0Lov7XSU6hmIGHSOYEOwRO1WDEGbztHbmPQoZtut0t4qJjmfOfM+mG65+pZ12WP9OPJKdYkC6R8vVNkGN/h0dXb515kM6O7lFpHq6mxD8CrwbJrJh1vahsDW+GXEPNTJQKZ4Wp/O0uBmU6aD5Mdi5EmZvld1jyzR5cWTNfumKv730Rqr+ntz8Ssjsh2h05er3folvfkHzAyPjIzaq/i4V/a/N9E7TtuHAzGOy1zEPo0DN5uaouknZrwhc/ZqzpDx/ZJwdZVO1AkHbb7Zbp5sLiYabX3xSBWyCcHFeM2i5LeMTaOQt/71lpJIlXy/Nmb87gMSGAsJIQVvcnkrA/LVenVx69a/3rpVeu+dotEr8eQDQ493MBrf7c5ufbr+yy0dkdBHVb58tPgxRXfUiI6Hf8btb27Q8M3zQ1OIU+qZk4q+e6vtGkuj/1r8y3b9GFL498PGtbUa7rEt72r5quEf9RjyyeaCT3+23X5DQEruU/7mL+vt4/PLv3KcwSRpu+UY8mHwD2/9OP/dna8+u/XFH7/7nmX41o71X9299Xy7Ptsob/MPD/UM/KtJN4T9H2pYbg71b92+bX9evDnUvxcjRN3rfqWJ3n6SWCY3h/oxdjp2ozNri9/TTvRDnKLXP9RXzRVcWD3WJ8mut1b1ZDmoEmw84lsJvS+xKrucamnrmX6xDFmi9Kn4kEvE3rx9j+N+DL5Y3xsXpXib0Xi+h+N+tLZP0jRzeIXE2EczKupSSp0rvTWlnnCNw9Pd98HkUS+HL0DsU25jUfA2Eds73DTqlWyX+5JzHaxSsitd6k2W8YFiUhf6INa3N/Bq8X0IMr6Ca4nrs3ENDINfiEPZmOrEGHxJbJ+DG0HH9cQC6P4bqbO9jW5qheCSvQk2j83i8NDnFFSxrBcVScr+CddJH4pP0yuhk947u6+IXRXoTE8krXD7j2qFhP0ruZPeFj+Nh1ZIyfsnSgfeZGdrd2MpIxO5T8Hw22QoE6P6nvlEa1JIvWsrCJd8a00/TRaeM87EvoxjYEKXZVT6l1I8Xn2h9H5SRjvc6EyPsYFeGRFVKaBzu6hwsTcTFqn3p+tTFYdwEZ+UKa/jbNRjyYLfUwWeHaZPezzDU6sXNz1gcf/s4/4TBC0Jvd83YqWLuff7VyR30fahzS7oHTtGZ/YEgV16UTsEvXAaTH2aVXQ+9NPKNfhV+dDL/pvqhtjn/Ru1QvajYQqrvcyeIAYDg9zGNPv4DvxL3g3v/fDJVm5dYOJ2vvnbWq2cRinNiVszCYp37FU3jLhKkj5IBuIfMuKvPKQ2l07fhS52qcsd7nEdTmjEZepwPyNkR6du2x2RaQgD01lLrGSN3twRVCZ1Niu9vNHgvaHH0Q2zv1jeaPDeix78tZn9MOPzp5+jSl61c009P2rqeOjVvPrt5t30xuo5UpZUJZ3pJXemd96nGV/Zmd7bgDOx6X2UIrHDQNBhvdku15k+m2L9tOn7zvQFb+PpmUBNTsFODAtMr3FG/UwP2eEC96v6vz7FlEY+QetRTcnE9Xa2F2w5J+6sE+wfQ95XWLiVVM06leW3vXEuT+w4HL61ac+/840cg0sHrRhn4/yjITqZKuC9C/0dmUI+GqxPE6ScErw1Po5Mf/qYWe09ztQj1cgKc3WhF/742MUeG/3pki71ccY5Go5wM+azVjTdY72ky3NuVN3/pd+fgmpkgB4n6n5+sMEfuoepdr3bH4CorgfT6eIJGee1fs72E/bUCtfvkaN9Y0K48Qkznihqhe9nB7Ha7DTz4xN+Pmyd+H52MOUk9pFIox8+2d7aXe5W569mi0/en/HG8YY33icmvOGNb3jjd3lIuuGNb3jj93TCfG3eWM2blPGd8cAfFG8sfXQu7WXJne0dbMjIisDwxZyTGVkFRMs+pj2vGhEtZ1vi+ETqfC+ppIl10AqbZWJNqQih5L00sgu9oWL8BoyWidlO7G3oQh/TTLRcutg7a8LULBXiiUMTYWPyrPnYm+D89PFEheQ4fTwiUM/OTwAGfaJYM3K1PCEhFxnHBXh8cmkSZfOKFBvy2AUPxLnESbLPN0r2ggJBL55w3sZp9N0/B+8sriccU6qYNF6dZA5T+6uziP73Za3I4wBzqfLkgIm1vo+Hr8S9U1N9JfVERMsg2HQRku2Iu3acwaYHCmmremM7xynqqHpCUy6e2x/cVNViDtjmzBNhEs3zhM09EvjxUmjGAlemWRlxbfzoTG1SPxpHxKkfLX2avwJj/e/XLK1uXGMPsllfC4vwT52mejzgZXHJBa9OpXjGplLwgc3FGRNiph5jWxaNHk8lRxs9OlWMbZMUXGZTiCYbQuwde8biu+ptzkZyNvGaY6wzfXAiJZXoktrPjm6zi2Xuk+Ugm7P1gdB8x+mw39Av1mYXghfvnHFhdHsdHU9/WhCBc+4ZiyrjlWm0fx+/2NvW9bm4KCmEnDEinzxXcYs9McUHXrHio4toQmMI6lg5dn6njq/JhdT7GE2MLmcHJox+r7jFnkCU8fZieduF3sYYnXM2SqjBu6b83L3zzUuVuHcaLeM1nWJjiPXKxRAf+l07xV4H9J27xGJEfsIltlXj/RRdsb0xMUacfHXV/EqXWHD5xiX2d3WJJUrKpEZ9JLV84xJ72kfg3bjEqieH1NOX7nFagQHR9QX2MbvE3vv+D9WbolOZ5Xq7/Z9z5OP3XDCp/lG/m9JeHVferkvsvd0fnqwvtM8fjL8l0q43MM1EoFVtM8HZ6nBZnYlvnGLfXhDPG7HhjdjwQxMbVse+ajqr+MnP1xYbntjZ3p7jK6AR5t0Hk1yJOOPmEEgo5G0KAebdR0fA7eRFekvw6Wi9hM7bmHtCSUdOkeTackKA4ITnnA0EpS4WMYgrRSxh4sWZ0GebnZcYHbHgi++zsdaYbH3M0uXMyTPF7DwCStdlxD8+55A8gZpsl6T0wdgSknhPbgkrJZheJEsxOb81w9Ajb8jfy7nWtdhWgQOdhvuD+fHe+d4Hn8TZktXfup33qzY6hhCDekgrZxST9DEk47Cutpwjbpxr62Y+OZV+BM61r+IIX+k0+aE61wZvbPbBe+9tsdcknG/qXHuLS3nlV44MfBw3q9v1WPzdvI2TtykFCckU/1aGY+Zamy0xPT8eR9K3PhYHjrW+N4Goph/peEzeCocetepgO6qa3olbbQi/j1vtyJ1ILjb1gZgg0ZNcxNsovXMppYxviZB4D3P/hK7HOs1GYk1fjBdbpERLEqTYFxuJRmgLWXAc26JP3mWfbUFX5KX3GaVliEK+Gcloj3KxOYWUSYvhU29SwhMg4x0kpvTRxZSTEU9mHDGmTwKrkoJLpCc1GW2q5OJSseLelbnUZxqu7a257u43G6JY/H8Pn2yPJRCvfuKfWi7xfPdkrdk7/rtIJV7lTAB7ycEhdmC9ehAQA6cjva1RlwHf2XjjJvA2EpPcyDRuZBofmkxDvZB+s0xj9AqYBBljxRu5CUzsgvPG+R47f3GBxFY5BzSYORljHRkYUyI+s1isT1LnYnA9cg/nXdKU3Cb7nsjDyBQ0NVuw+AZKhMPwpCPDYTP5YLMNRpPb4SpcXLHRJR+j7bIrvRpaOReit0FlF1l8LORvS9F30ePBmRz29iaL7WLKtk+GM30qyb21UGET46j0+N2HQ5km4oZvWzz+5uH/XT/a/fHTFrZrihemjp/vlG875tlutEazSG5qwH5c/vCiRbD332iN3nVIjxsO64bD+tA4LDbRD5fDctn5PngJUdOC5pRMn3LGbS4X6WLCtFdicd6LF9e5aDKsDcyTC4akxj72kVAwyVkhNavPuTculWKxEyfJZ/a9NcG4BP/mu2IIz5KDM/hnOttlk/ucUpRsQkZAEwshMIIj7qyllSgWBzZnxBlHfmDiy/Q4gMZig3HVkX3ak9khP45IIjcslsYKr5ZAv8hivfW4dp/e+f7Ord32xe7prSYMGkNr/NMzWP/txF8cnd7TxvCKNLKvFwb9JlL9+0gmcBOpXind7JTUMj58OHEHX9uXz0YMTCYOK1d3gtE0Z17USPUhxd5JtIVUWhjw6ZOErishOPJ5heS9c56heMfhMKa9P+Fi4YlHIJjkC5IqOCeXSmdD7AQuihTv1nUpd95//BKl39LljyVCwcOHt642Ty5W5/+Qr+C8+s9k7nvDV9xkwFGB3BSDv6Wy+eWIxjdimxuxzXvizl+br8hk5ZyJboIzZFpscYLydIccNwG/Tbw2iSlk8darXMVi6UNSV8ZkQiYtynTLYv4ZTSCRSo6p1Iyn74nhiD73uUQNlkSI267YPmdpuVYsNjTEpR1jJVlUb6SBmaJWRdfhzz15PucSiS4whXY1cDBIbKY4pkVDysY4c9r2BOTKkWheKOvqyHIRFSykWdxQApr2Jdi9m/pHzwBZeNAxhhRW1s6U2JsWvyzA7GmE1tG33LrOF+mnmASpkIUo9mEf8YrgUiThmAaSYLPzoK/WdYHwCNMDKs/r7T7Iq5p/k+9nuqST0ot0MfRm9mK6cWvfHbqw32R8PuHzrC5+JzIE/4Jb+4k80IcJn70nQuA7cmt3uc/q+I5nfFGnxN0+ibLzvVBPtJBr6Z4Tsb4JAVGIaJLYBX7YzZ3afTIfulM72xFe7TaUbFyg/6Pb+ps7tTvpo4ZjRE2QNLzB6LX+E6miNebGmK1Z84WPtzWZtjj8VSzhDpIY3FJnXu0xNed08mh7QtK/plf79XTO7zbV83VA36dXO177P7yseTyZ3DHRM6P5K/3YVeayefTTGEZ08+inn6sU/8vLJ63yy8sn+7q/1bv3Vmf3VhuS/cD1391eXl09XW22R8V7l00lMHf+hec6yIeLiObybPN4sz5bLB+vzq/W3eLxv203Zy8Pi7PEo3ceXt29vNzOnv/sbLNbPTxfL8YcxZ9vXq7P7lw9Xz/afbvabS7335pnLz7yLPqlDH+HZjufvVw/+nUZAg/f36d2rtm2SDKyWOLHruUbP/bfw4+dhXXgx65Lao68yBMPkPetS7M++eY+iqvXMR63vy6d87PVbrvRLIIjtl3TkP2ukiynOva367l+I8r6WERZH17a7TeTZ/3DxADvUZl5oyn759GUHUu0bMilb1IroiGOsilii6U+ljE7DNmT9xIt632UPktxxmVzmOzZxhis9Na6VKI1bIrvWKKllB//Fey6vWsJbLKNSDqcN020EchwFPqyz3kUCGXpigbV1CsSUtOkuI8QSAT6nGa5mKhwpJOdLmRXJFaaKvQZQeqND31M/MdTECtkjsaBu4lRcIM3SGP2YhWNY2+C5YYxhVj4xiAg0x828aP9z4Xe2JYd2diYOtcH2WfnEWAVcQTGx0p9+lJv+E5Z3/bTpwxx0qdAoiXnQjIsvl5HhrEyvQ1lFktdISTHDB/UMP7TfxqYUxOIAWcRHSOh0oufhjkRlz/iQzxeCnOaJbQivKPpjSQ6ZYyL9eOthfl/hNwB+wCTkWxaeTrioWBNvZgyyRENyQWkzDIl6QROTbUpoINjT+fNmb7EcQFVmafpk5tJ4EwmhCnZh9ql0s4S4n7OdaYMZnF00et/QJlp/Ez72ac0y/mkQlJvxlxaXMhTvfHjrKlDRB9lj53CdImkOpqtKccQaIttKrXhsfU+JJkyKCAeBBFmeKCI4BLK9PFynSVt976GRRKsnS0SXR22wsHgzpqseDQ17+OI41wMjqRD8aQhsdt+OHX+EM6M4DBStjcyRZWtc5pkQoIRv8ZlkipWj8tkgsXl+VLVpU8EOIbU0iNrZsmwEJyqneR08YKLcSZbpcZYzXRRl3lbRHVl9i6oGWa7tL3akEl+tqC0PyKaJ8PH/bBXhLBhCoVLuFJDar5KVxpuh4px4yyA587MMoXVpV4ws9hffCckxZxURhG1XgyHswBbykw2rVhu9wn3NLUCvVewD2jSSElsmucB0zlJU9DesS3jReF3WUmcLi9b5vPAiynYvVRcVw62qLQ/dp9VoH23Zj7JSnBCpUjO6Xy1/3jto8mz0a0zEeM89RhQIemnLbF1qdW2xm5KkRmu1gbrBmCL9sl57bv12nztr+Q5HvFSkLkYXxeYs3VqmeGGxWQvqR8N2lJbg42ySt22bN01guMRH5paSC9dgJJ16EYCMqf+kudzpk2MLdU9SCoOtnfrqxUHxOvXnZeZOkK7LEH/syeKvdg63W1MRHfSqReVsLXRrk3gK2Z65qqtJt0LDPkb9+SyUQZTl4dSAXF1i3OiI9oGqP3HZO3gYckqKZgCVLdLqcMsmV4dmjZyirQm5+N1NCKAiE5Z5SMmQqFchUEQCipVUJSWn0IoncrxqYMeVyowLsNGE9qecACf7uv75uuT0emCD6kOTx35Mimt9AK+9v02UQ3aow7VYRgrK1mqqCjR6wg0fiNPeXD2328cTgO8UfA6RON/4EsqpZnW3shjRctMCpF5DIEB5pGmFf4Rmdpoa93Y2LhPzD7emqo9qigvIaLETL4li+QCrV1LsVovakZycYLk+gpgI3at021sG7VolRAN6b2fhemuTtxvwUT+SPL5YcVUstJnZVichovgUPMLMZVsgfWCGMRQctCIy/L1p3dRR7Q9UKIpzgaNdzveiqXyRCn5HPWYc9+2txq7IsUk5zU+7IPplqtvWZe96K1j+eAfLza7GtLou6ers8u/ztza7l5uL9bbb1dnmxejsfzzljb14Za2yUy8eoYkmFAPV5RfN3JO9VBZvxxD52iAxNeJqYTaYK8FNVg3YM3266TbR+j0W+IG1RCbm4snc0nlK8MGwQKpwuMDih2Uiws2JjHe+/y2gynVg//HFUwpWKtxqnMOMb29aEovb/3rrVsYSYRfs0BUnrdfH2Ox2k9Tul/VVw/qf2aRU1hKrzbwRKH15ebZZvf56tHucluN5lCsfbPdPNlcTEZ072o0fmY0bpU++f/mozFGl1L0cH0yGobn1SRkJK4fSOSxt79YxgFRDHGo5j9iDJn0pyfibWnQrXcfeaukdxR567tHq/P1YyUfLRgiQUl/XixBitAXK8nnFNMkwv1ZM0okIe9gSsFHtVSp9oq8lkga0GfyTfqsqog72+3q58UrIwhNIlhJkWTeenwTkzrJ3nMc1iuH2AXn+n2GypA6SWWegzx2scyt0chcXOZJycnhaZsNTbu62M+PwkbIgzRL6qJJPJFo9qHz/V7KydVZybkP43mI0KaaG2h/KccOlw5jziEBLr9UUYH+FYlpn/pGcieS+pnEzHeuT/NcjRwdVACjR456klABb/uoDa5l+9CrdMH5McF8KSnELvfVcLNdyJnquUulC/qnQdi3L2fp5+de17mS+0n8lXMXfJ/3p25BgqfShsCRSaWG9STeJNHTf3w6SPfeEXh2L9+yNVH9TBSXu9iHmbSjSi5opMpWVAjT1z5MjTg7md1yYfbY+5mYC6cZj5B7vCKZfvZST16RPnJO9vypR/lpzOmaZdrH19UNxx0kEyqdj33e94xM22aOd0xz2Uuw66m6tobwpjapU6OiDBoXjh5d7meCKRUa2v4w+VIpvZuy2rIAUj+TlSs+zUTQo1ztoK3ZiILgkzCat0OPl/feDKp+wfqDGiLakNVrvAgO1qsgqaKyypn0ULzvHvnJ9gupdLiRz8qFhXG4TuMcCRmLkFogXL2IS2YR9u4vVdCMQHFxjK8Lq8nzVNkwP85fb9b2BBEkI/DsYqDm41Q6T2KqfVtk7VURoGJulZSMg96J6S351w5FLdgwW/K2+c4fJq8q2FvP5WxkLZuLEbVjVaoxDjOUKHZF0bd3yP6dwqQ5irtymNMrg+IzkTyLBhBbCmaQtDOQ0hmdUMqi2YnJLF//eKy6U+ysnaMDU1UJhooq69eALvYz7Vn9YAJX5mK4Wq1S3Ew7DiIxl2UCR+rNXP6shuMNeMbBdaFP1Kqk7PpyUBKj0qDEGBHxURPb7S/W7+F8N2pRJ5WumP36bzBDuHQBREZSlzs/kgXZD4m+yFHWNSS2E2Ugq/E1fZeK9FQMVsFXAgY1QY7kbFf224JezAzLT7cqRbE53uuYId7SYdA/KqXa/0qhc6zB+YICRLpkmXcVgx9PnSGtW6NwVYpGIyqHs6SpJu3fdSRg5GzqAFWxui4ilXO35uo+Wz/p2JCESVbhqifX4UFGZCDQ7+hwab8UDAVoD5Wu0oRKjsCi7JqoB8dJ1OaVlOm7lWbratM/wKWAGNIk7tMbcpF9e5oe7Ts9q/3Rlynuv6UQRfIOWtsFoq9PdKURuTqj9apr6wQxIZlfr6oFRUIdgOvTE3Vd+M7nOYURuTHgvzHg3351ebauMov7JKCzvU+klYsxeBs8QjIs+CHKzppEQN3gAgIKZKq3EbVaHg9R8wnwOO5QfZBirHG+ZFfz1R1a8OcEAZNgJTqpqdPG1HM/LZY5+t5F7qZMID1aazb4i6XzuS/ocQIR87Pmwjswsrcp+T5JlJC9k4j8ZW7BzxGkFO4GIgTO7tKfSCqElHJOMVvNendowR/n9vuowd4kLZ2PqbfKdYrR7H1/mnka2NgH5TasAjneYo6Obu0Tx2GizxSWJC6TBDUUFcnOLPhPzfB4G5m3MTZHDD2ykeiKpibcN4CZvvRklrMSTLImCxMwWunr3WuYMt5dLG+rAU6IPjibQ0qqWD+w4kcfo1dK0WcbFAH3qetO3b7/w/MmeXdVYOhcCDEHndrpnh/5GmtdEMT8P45vSUuHRz43FQt/0e5MZwzvU3VjvfPwauZycAqaQ1v+Vzzx6qbvPLx6des/vHy2uagzrMEnYzSJFMDVjJ9bgkse2JslZ9F8Es3Y/7b05cRLZLSTPpjgXSEmt8lRPRSfrnerxdIh5rr3dNN+HVjEc+vrF88errcLtX2+8QT4XTPaVfnt3hNAyzeeAO/TEwDSerzE2uZ5vCjb5nltRbbN83g5HmsJ360XAXLe13Ah+O75+eZirR4nP1xs/vPFmsCoUIVPNrurhWYNbOrI2ZMjhp5fPvnDo6eb//W/7P/8H/zerV5oirjXVpxA3T/fnI9xP2mX4ne7n/GF4e5Xq+1P6+1MgVorrj1RdSy88en6HKp3m981UuepKa353k7R06+fN0L5v6926+d7p6EvNle7yyfb1TNEqP+x+3lPKb6QzxntU+MXrFLYqtCtI6gC+LbEP8gB/PrR+vz8arEUzWL4dr0x6Pw/Tah+pva3YNhNgkH0pauHN8H4rxZLdFl/7xY3Hhr/RB4aVkmD6uLqZqo/XzvD4Bh8/y1H4zcWPJvMjSDFH3JsVUd65/M15yE2oufKD+AignW8JOKr8YPAI/ywptWgKa01jngn3CL5nv5QIRr5zhIyN2+caTVuvOW8ipY95tXtmWLr215CrUEgq6/7GNqtgsCv9yZgEKY/fPtySA3AaEhT0HsTXa4PxyRq3WcSOZS4lbwqubzRbEaElMEJXn+EaglocvH1M4XYIugfsoZV8WI4aAs/Etkeey9iCSbCj6RwiiWqDRW2QSXO1WETp/4i8L71+xLauyHV/kt0SaGUhCkmj6akugov5I/TmiKlPlNK0Hmxpj1jJXgdGvUq4mHrrKs1nhwP6vNhSNqE8wceFvyIxeqXbSpSazJB9FR3VufZifHaKWfbYDlXv+K8J0lU712IKIV772KugLosddRcqUiik8Aj3po6qZ6UmFR4TTXVex9tro8Q4lhvlTqcweQKXLBOddRqr6k9CSFWGEJyqJWJfoxGVXofsTvUHzbU+Y/eqTCYfBUVZUnmqU0hzdfGk7SRTa6ousinUCqSoKDWJrIxqkbx2WrIwd5nTyREfsTaVM6pDkkxJeksYB2pg1RCCApgSbkwskHtKv/975CPf9uunj+tuxas3HMkD0rt7l2e/7xYHlsBVpaXe688YigXdVoe8gCxhjcSfTCC10FESvKJSrhstjYHEoaB1tiHqrWnmFxMSNl6XzzGnhxnSKOaCndKcVbKNTu0XwUmobZJmTqJbRqYFu+n5B0ZyCYwHcsea2ExMRbqK5heCpG1M3OlUsN3AGcV6l6D0yEX9CI5aoCNOpw+lBRNCFlTtCGwq3AmjrEmeZK9aSwPhZNFlyUhPiyWfeXoRPTrxtPFQNSP4/HUyOSlOAnVyFanPZB21lpqHcbRI5wl2CjGGWuiVzvYBmcIAVFmANLr9oe/Dk5ryql5D8GzHWXjVYjX4MySnfigeWK1XsfTSohQj4isOIC1FU5cXbwvFqe2/GbDaXOJiP6PhzMFXS4p+MzRuIKZnIhI8C5am2VaRTZadtZorTdJR1/BDCHlFLxSeevdG8IZ8qllJEa8I3GgC4ImosJZRCz45gjKo/JkHU4n3knxeO8Yo9UKZzSE9jE2xRKdf9PxtFWVcjSeunO6IlbKHj1x7ArOOZuMJ91gw06XSEeIi5WLoaiQAmSO3tiUJYj1JfvwhsMpuCZcm3YJErOkQD7miSqJ9ZiFiStGqmZGR9MbCSaIJJ+itzMwURZplsXonbwZmPjvnVhEzG00NlljlEjqrItL5E22xkRrHPBUMJ2kwHSkVIKbaHxkjwg2FrJVmzdc65KNkxOjWXzySYKQN2pEToZYiqbQNkEDGlUwo7XwlzljcYA8vk06eaJ9tNEmF950NCPpq65NOkm4i8noozVJdR3NGIOLIQUjYt1E4X2xPkj0ln0rQlArmM6EYEuKNnhv3xRMX8wJ3CTJOBquYoKG6KpgZutyCnh9ZKueHDqaQdOSZ3GssD1FiibCcnEgMFbecAmJO7mx28TcmmJSnC2hUoQYVSmWEtV7voIZjLfRRkdQSFXC6GiGhFOS8xbi7K5nFf5V+5BIyidw00EkozESrAZr1dG0EpMmaQULPUhYwUzJivHWiMCOjpMeLIAa56y1Pr0hfRejOuVjuqlMuQ8Jr/RxBaHMyzEIfmI+TQs9GsmBZCTBZ4JZjptlLBYsTtjPufhm9AhV2QnUdOiNS/KSgi6VOpjBO5fFmYh74LRZgoHB+OxLMVbjBemcc+60DlbeOFPecDANyaivL3RvXbE5EIktMD4VzOTFZktOFsPBss159DaRF4ZVFyzKfwXT2eiMt8nEkuKbssYmhFNkk3XuacdV3qyCWbwPsHI2Bb/f0mPIRrzNJOD2GrFCwbS4Wjmb4QNcfMMtHT/SE2QzSCiRaSfK6ziaToIUSxBVToJ73IzKoNjoWXLjYEow2SW2yxzkTRk5Y/0pqhlizCFY9kU169bBdDBwuDLGjLn2NOds5FZcKE7sfg8yATa6EGEvGf+mjJwxJxm5CIVxhCkouqQrnCG5kmHRrXAeHnET/i2H4EMJ+MYrbt4m5oZNJWTH4QMu+03WuvSct04QzpgQrqSQc9aTWAU05ZCjg+EEgD2gpKzE/09cVov7T2tQQ2+Kz8laZVPfEMxmu3PEcOp53zsnmcN9W+uuFOw7kBNh67IH07gSgs2RM+l+PK3TkBVRDOPwRjRJSL95inSy0UggHEpUEyMdTm8N5CuwLHza42coDsLpA4ZKezidzdGYEp1zKdg3HU9/kuXMgOiMWBfVNqfCyXnHRMvWDoOnu2X0rmDs4l1JeVxFt7GYzDmk4J11xps32oikLy6TG/141gsk1YgtmDONs+6DSbFYl8SmvD+sR8cpXgxBV5yaPVXk9BYQc/bWBV/eiJcDzpOccfFCnvXMimV662hGE01CcmMcnuDjeFofCT4CESYow7jafXKG3GEc2oN7IxqPibYN4NzRKioJ87eEtdCexvtkbDYxZyTFeyFNFDLFRutiMkj8TsCZXHhTonTN2KDa8gg5aNlfagjaOp6ppCxweMUm5S4VP41JeaKaGH8R+NUJLKjLPuI39fhPq/NP1k+qNuJPq/PPLs6amn/95O7lWPjs4qwV3npQw+eoG17DHqHKF1+hTRfjIHFNnT49OmrTp4q3ZX9gVLJ5YIKAIHm0SmhGCGwOzSph87dmtvD188vNBQYUYM39xfLPty179IxeMikPmjZjLunSBFcvLh7hFP7rHLAn9zjNKjwzV1B3ucV1ddK1mF573dBbn/+H5825/TVwgEgE+Gy+Cg3Y/d7SFKvwuZrMEi46FjWYRX8ymswiN+Rkha1jNWLFZNZiFOjFpURYqmph20xmM0xujvA3kLYpDgEyyBCUX49RifYUhwBOzJuIGAsOjrfGYAPJZ89CTuKd0zPQPg5Bi17Azo9qZoxL2tyNu8VvjUMwTlVbZzs1VhldZS+eYJrD15mifVvqZFv9aC2wUJ6Z7VCc5alS7X/LwPzl5cWT9dVusVRbmQ8hDsGrzHK0j68xAOPotP6PxffmZx44e3qLcDGZUA4dZ28cBm4cBo4cBqR3KcYixmd0p+oAIIul70tCoAqXUypbB/GTPtjo8SqOHo0RpNIuloRhjAZVZkacDRE7sOi/LSnG3mIwHvguVHx0CFCLd2NCb4XbWQgnuL+9WKJJ6YNqLaOq1tUdYDRpJ6a9+NhblG0uxKBB0ef+AsUEUmzEXHm5+rqatGPP7jJqWewCJGeVmB76C5BOoLqVQYyVA/hhb+TP8h6Li6WV3Dz0UK04mvrTwTi8E38B6a3Y6Cy8YXCRwRv9Aao5//H8jncXy9jnGIJ3yeUSXGzW/OPQ7tRbIPTZE00u2VKSitxHd4CfTqLJeFfPxoQG1PCDCHtCzcYws793vlnuiylkzTsI+X/K9H7vKlD5qRhjLkYU5d6tq8B1UN+DowDeAIIE2ibjyCKo4f7xEwjESnM+SJJQVOnZ3ATw3BHvrfcqymKFtgwBoXcBkUZx6kLPKv2VXgJwRjdeAr+rl0CNebP3EtDyjZfAaQtskP8tHRakSno5L5xYlOyY11dk2zCPl2PbMI/X4vt1EVidv4dD+Ze3Hr3Y/kWt3t/eNLDlHhzLOWb/Q8+A6VCuHlH3fyEUiun2loZz4yOskQ7MfKg4sKehYn6c14oDCxGeOLDFwJLuwOiBCuzyJusCKg70+Gp8N9eYU3Ggm8ZC7kALTMWBvpWKA80m9nEHOkQqMKEbtXWUD/Riai4310Cpldxc14O924FWhYoDBQaGcAeaAioOJPKCYLC4vewbM7YDITMVB9JcwVNwLjil4kBCScWBKBAzt7nMTYUoD/4xjoTeEYC9kCLPYjUReiewmTlmUEfjkxuCEsZMbIDAA4ZY0pCT7IvVCkOY3Ggszns4lOeCUSEGfV7LKflExMyQc/E8EL2JnPyzcUR+6W2IWY3MUEABhcXkxCWEhcXSqvVijSuSkTNQRp3mSvCSjdcKi0GKI8RziIBlxWeXi1Z69XM3wdInr+YgXcCcA7F4NOTT4gnJPrtAhMSQJFGRvLGC+V7BRTgQFzOzRAQ5pT6B0NUK/I3UJ1yRIjaiUfH6io0ZdRc2f9qIoBnzFlPaDKAYOsTsYzI5GSpMTsF5413AXJacdziGpiicXLQionPyqcSYglaEkmxKqgbSb5pgA70rMeeoT2DViVOqoPxT/HjLgrojK7pf8itSdu9VgtDTNDeq+AgiPPOeovjj5mz3tIoqf8EZq25+v4Hq3rgTna/nkjAYt5kobEzZrptadSc6wV1Uf7UTjP2NM9EH6Ex0gg9UJ+0Tx6/3P383uWBetQQROenJGv6yrTh+PdyginmP6XpeO7uxipaJFPhb/IxOYOPkcnQdg7++3D6r8U1rW3cvL3aXL7ZtXL578QwHJx2jvRrpY3IxUveGqf+TD9bq5ebZi2f4E4uaCew799Z1ZP/ddaTHZ6dRRXrMMN+oSH8PFam6T0XikhRbkpr2VJGHS9jMSCLkjNrwVJmHx1o8lRy8x4dpryTAoS4Fk6o6c9KPeiz4nU3Y2lk1RN7rR1twRVd8zjVl7KgfDQW7yGRTTBlPM5ppgV/+6fSj2O9/COrR2fHiY9GI+pjUaK8lAbjRiKI/MD3L8tvNkxuN6LFG9BU50E9UQ+xsbrF+RS221HjkdPWhSvRdpkD3JfelZFdSSUhoZirPn4gqnXubsisYk+aDu4vl7dLjVlFKTs6h8KsKzlEliMrvbaZAfycK0XebAj0ERgghFznx1GJnVHn+Yn71qhAVMc5HaJIr1cVlpg69nlj8tSOnTUEjp+ho71Ydeh3U96AOnek/9eeJLOjU/0qtpr5ykwX9ogYq+l2yoEc1MdtrNbV8o9V8n1pNVgE8Pf/VfU0r2Mn48X6Vk287C/rdy+32D7vL3epcQ5SNiAbSzfSSv7eMXF133m7IrXu7PzxZX2in/2kib02COsVUFa0qir5/Sd2D1cvNdTPtO1S+YukCfFPc/Kwvq0HNDqvaD3KuVMQIfM3Wt69BVLXqddQaH85c3Qi+ddJO6J6mOaqzpfFkqftwBd+VUv5m0fcozb0m730tKfd/Lf5ltd0ulovPzy9Xu0iU1hqoyvu3FGLrSBP8geVmJJ864hQbMUhQS9163K3x72MIxG1oIj9LLgJ92kqswcgnmd9NbsbZYjzpE/HB5GY85p8+prSMrmXxkJR8CNfDP/2WXJWfXey2m/XVrX+9Jd7funXr40rLGEuwpgRHhI8TYTN+y4B8tV5d1MR74vtA8L6PRjr8zoZD0+7dDMdiOSUi1MyMt2LvPq5Ule9gvRykZvzoR2RyPv39cjNKzu87OaOejTRR40Fw15+rhcQ/SLioLKwGNCWVliMwaMihuJCT5j+z0RLNSn3Ag6YtskmM1WBhCUcQ01uH84t3Sb2qyLvkjXPGlSTWORN9d5v8RSkTPSA7gz62u+16S1QL4kzgIZ7WtzXxW3Q+SbZGHbsseZBsxuDRGSLfEeHU9EiWjWRbSOmiX88m42gSjMazPNETTY+VkHfjDkUcQT7tBFNWH4vxLmuKLOsIvUnQKZOCJpJyxVsnUggf5aOrfbFE4EkRqHSYjJHirPMhEajNC82ZiJ7YOgxaiYwB4CkVGw2RGI3RLPPWmpDwQ8ukAbd8PeGnpBGaYjaW4bw+CcddEQVLMFVlmHBf0/xXoRCqI6bifGAiDEaqLuqQFmIj8Z6atBJ4xNOdqLBjXVtcDDbEUILm9TLWY74aovfiomP8iBBrnYsu5SgiIIOJkViJ4jAmJw4tc5/xqopiigNPGPhrSHQ8Fbx21J/aSXKsBB+JSUpgVVLcMznEFo0lWW2QjPUmG5uSBEcKMeK6RJOCEM/HEN0JCYLLPqRYRB1UAZ64nVgASAneKDoWjQFiYyImR9L5sVGIrIUlbyBwxKkFcB2Nrk0FHTzqS+0gUXEIk2G9J6CtAg+D5H3OOWlWNiMme0N6b0e4On3IEcjVJC+ZgIM6+4raYGQgHI+mDUtYVot1Pov1EtrsJC8m8G7UfGPHy5cgI4eL4DoenZyH475oB22xxeFNVryJtTvBZOLalpKx9G6TI04kB8lFkmgGNKzGSWheIqbhtEgQYBtdiESz1Y87KZr/jOBhxLM9RXkgRocL+NQioO4EGs1nwpzsT4XDiNqjk/mK1FCmJzpvtj47gr9ila4NFE9UQyvZWWzEmWy1pg9SvMHERXtuI4hKgGbvwfBX0M5TROjUQj5aDLR6AqGOJ4VYM0e9ql0NtniTJFnCFdVVE9EtErQvl5CULEgkSRYOA7nEpGuSFEtecoiRkJ6Ksdep/wkqeooQXV/LpxaEu4ZUp2bluEO1l07AzpBzFGeYKXHGFONSVov8ovPpbYHUBU+AIaWkJ3ava5vACTJ6nRJdX8zHC+I6QtWEg0eTwTwedabiMRMhwZeSsovaH7aZYCEohPkCi4+2XioON7ATOwDZCo/IqG4dh5To+nK+vhgUR65j0omZOO5MnUONDIRDSilWdCs/YhxObL+Qi6P96zYZWQ83gfrUAR1lnI+oEZTncCmz+x4tA6bxGI/4/NFc8Pmj/rQ0jFgmfvdou3lO6BfVLX2xefL0HAubu5cXF+tHu30+mc83L9dnVab+eHV+tf77/wN4nztC').then(json => {\n",
       "   const obj = Core.parse(json);\n",
       "   Core.draw('root_plot_1779221491743', obj, '');\n",
       "});\n",
       "\n",
       "      }\n",
       "      const servers = ['/static/', 'https://root.cern/js/7.11.0/', 'https://jsroot.gsi.de/7.11.0/'],\n",
       "            path = 'build/jsroot';\n",
       "      if (typeof JSROOT !== 'undefined')\n",
       "         execCode(JSROOT);\n",
       "      else if (typeof requirejs !== 'undefined') {\n",
       "         servers.forEach((s,i) => { servers[i] = s + path; });\n",
       "         requirejs.config({ paths: { 'jsroot' : servers } })(['jsroot'],  execCode);\n",
       "      } else {\n",
       "         const config = document.getElementById('jupyter-config-data');\n",
       "         if (config)\n",
       "            servers[0] = (JSON.parse(config.innerHTML || '{}')?.baseUrl || '/') + 'static/';\n",
       "         else\n",
       "            servers.shift();\n",
       "         function loadJsroot() {\n",
       "            return !servers.length ? 0 : import(servers.shift() + path + '.js').catch(loadJsroot).then(() => execCode(JSROOT));\n",
       "         }\n",
       "         loadJsroot();\n",
       "      }\n",
       "   }\n",
       "   process_root_plot_1779221491743();\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TCanvas::Print>: ps file testUnfold3.ps has been created\n"
     ]
    }
   ],
   "source": [
    " output.cd(6);\n",
    " histCorr->Draw(\"BOX\");\n",
    "\n",
    " output.SaveAs(\"testUnfold3.ps\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac136884",
   "metadata": {},
   "source": [
    "============================================================\n",
    "step 8: compare results to the so-called bin-by-bin \"correction\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "578f128b",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:31.757330Z",
     "iopub.status.busy": "2026-05-19T20:11:31.757198Z",
     "iopub.status.idle": "2026-05-19T20:11:31.961373Z",
     "shell.execute_reply": "2026-05-19T20:11:31.960578Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bin truth   result   (stat)   (bgr)    (sys)\n",
      "===+=====+=========+========+========+=======\n",
      "  1  4002   4261.3 +/-137.0 +/-  8.9 +/-473.2 (unfolding)\n",
      "            4190.1 +/-128.7 +/-104.7 +/-266.3 (bin-by-bin)\n",
      "\n",
      "  2  1816   1775.1 +/- 82.0 +/-  7.3 +/- 25.2 (unfolding)\n",
      "            1714.1 +/- 61.0 +/- 17.0 +/- 70.8 (bin-by-bin)\n",
      "\n",
      "  3  1011   1020.9 +/- 60.1 +/-  6.7 +/- 24.8 (unfolding)\n",
      "             959.5 +/- 40.9 +/-  9.2 +/- 43.3 (bin-by-bin)\n",
      "\n",
      "  4   593    526.9 +/- 45.2 +/-  6.1 +/- 15.0 (unfolding)\n",
      "             506.2 +/- 27.9 +/-  6.6 +/- 45.2 (bin-by-bin)\n",
      "\n",
      "  5   379    371.6 +/- 38.6 +/-  6.5 +/- 19.8 (unfolding)\n",
      "             351.8 +/- 22.6 +/-  6.0 +/- 29.2 (bin-by-bin)\n",
      "\n",
      "  6   256    341.5 +/- 34.4 +/-  6.0 +/-  9.2 (unfolding)\n",
      "             282.5 +/- 19.2 +/-  5.2 +/- 46.6 (bin-by-bin)\n",
      "\n",
      "  7   193    182.0 +/- 30.1 +/-  5.7 +/-  7.5 (unfolding)\n",
      "             177.7 +/- 15.8 +/-  4.9 +/- 22.5 (bin-by-bin)\n",
      "\n",
      "  8   137    147.1 +/- 28.4 +/-  5.6 +/-  7.0 (unfolding)\n",
      "             133.6 +/- 13.9 +/-  4.5 +/- 20.5 (bin-by-bin)\n",
      "\n",
      "  9   123    109.6 +/- 26.3 +/-  5.4 +/- 10.4 (unfolding)\n",
      "              98.6 +/- 12.1 +/-  4.1 +/- 20.7 (bin-by-bin)\n",
      "\n",
      " 10    78    100.7 +/- 24.7 +/-  6.1 +/-  8.4 (unfolding)\n",
      "              89.7 +/- 13.2 +/-  5.0 +/- 17.1 (bin-by-bin)\n",
      "\n"
     ]
    }
   ],
   "source": [
    " std::cout<<\"bin truth   result   (stat)   (bgr)    (sys)\\n\";\n",
    " std::cout<<\"===+=====+=========+========+========+=======\\n\";\n",
    " for(Int_t i=1;i<=nGen;i++) {\n",
    "    // data contribution in this bin\n",
    "    Double_t data=histBbbInput->GetBinContent(i);\n",
    "    Double_t errData=histBbbInput->GetBinError(i);\n",
    "\n",
    "    // subtract background contributions\n",
    "    Double_t data_bgr=data\n",
    "       - scale_bgr*(histBbbBgr1->GetBinContent(i)\n",
    "                    + histBbbBgr2->GetBinContent(i)\n",
    "                    + histBbbBgrPt->GetBinContent(i));\n",
    "    Double_t errData_bgr=TMath::Sqrt\n",
    "       (errData*errData+\n",
    "        TMath::Power(dscale_bgr*histBbbBgr1->GetBinContent(i),2)+\n",
    "        TMath::Power(scale_bgr*histBbbBgr1->GetBinError(i),2)+\n",
    "        TMath::Power(dscale_bgr*histBbbBgr2->GetBinContent(i),2)+\n",
    "        TMath::Power(scale_bgr*histBbbBgr2->GetBinError(i),2)+\n",
    "        TMath::Power(dscale_bgr*histBbbBgrPt->GetBinContent(i),2)+\n",
    "        TMath::Power(scale_bgr*histBbbBgrPt->GetBinError(i),2));\n",
    "    // \"correct\" the data, using the Monte Carlo and neglecting off-diagonals\n",
    "    Double_t fCorr=(histBbbSignalGen->GetBinContent(i)/\n",
    "                    histBbbSignalRec->GetBinContent(i));\n",
    "    Double_t data_bbb= data_bgr *fCorr;\n",
    "    // stat only error\n",
    "    Double_t errData_stat_bbb = errData*fCorr;\n",
    "    // stat plus background subtraction\n",
    "    Double_t errData_statbgr_bbb = errData_bgr*fCorr;\n",
    "\n",
    "    // estimate systematic error by repeating the exercise\n",
    "    // using the MC with systematic shifts applied\n",
    "    Double_t fCorr_sys=(histBbbSignalGenSys->GetBinContent(i)/\n",
    "                        histBbbSignalRecSys->GetBinContent(i));\n",
    "    Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;\n",
    "\n",
    "    // add systematic shift quadratically and get total error\n",
    "    Double_t errData_total_bbb=\n",
    "       TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb\n",
    "                   +shift_sys*shift_sys);\n",
    "\n",
    "    // get results from real unfolding\n",
    "    Double_t data_unfold= histUnfoldStat->GetBinContent(i);\n",
    "    Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i);\n",
    "    Double_t errData_statbgr_unfold=histUnfoldOutput->GetBinError(i);\n",
    "    Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);\n",
    "\n",
    "    // compare\n",
    "\n",
    "    std::cout<<TString::Format\n",
    "       (\"%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)\",\n",
    "        i,histDataTruth->GetBinContent(i),data_unfold,\n",
    "        errData_stat_unfold,TMath::Sqrt(errData_statbgr_unfold*\n",
    "                                        errData_statbgr_unfold-\n",
    "                                        errData_stat_unfold*\n",
    "                                        errData_stat_unfold),\n",
    "        TMath::Sqrt(errData_total_unfold*\n",
    "                    errData_total_unfold-\n",
    "                    errData_statbgr_unfold*\n",
    "                    errData_statbgr_unfold))<<\"\\n\";\n",
    "    std::cout<<TString::Format\n",
    "       (\"          %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)\",\n",
    "        data_bbb,errData_stat_bbb,TMath::Sqrt(errData_statbgr_bbb*\n",
    "                                              errData_statbgr_bbb-\n",
    "                                              errData_stat_bbb*\n",
    "                                              errData_stat_bbb),\n",
    "        TMath::Sqrt(errData_total_bbb*\n",
    "                    errData_total_bbb-\n",
    "                    errData_statbgr_bbb*\n",
    "                    errData_statbgr_bbb))\n",
    "             <<\"\\n\\n\";\n",
    " }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a635d645",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "178591ac",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:11:31.962935Z",
     "iopub.status.busy": "2026-05-19T20:11:31.962812Z",
     "iopub.status.idle": "2026-05-19T20:11:32.167465Z",
     "shell.execute_reply": "2026-05-19T20:11:32.166772Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "<div id=\"root_plot_1779221492071\" style=\"width: 700px; height: 500px; position: relative\">\n",
       "</div>\n",
       "\n",
       "</div>\n",
       "<script>\n",
       "   function process_root_plot_1779221492071() {\n",
       "      function execCode(Core) {\n",
       "         Core.settings.HandleKeys = false;\n",
       "         \n",
       "Core.unzipJSON(88846,'WkwIy0wADlsBeAHtvWmTHDeSoP1XaGmzZr37BsPguJFh84GipFG/q4Omo0Vu29hYkpUkc1Ws4mQlu6ke6/++9jgQkZFZSTUlHiJ7KrpZSiAOOACHw+Hnfy3+Y/fz8/XF6tl6sVx8f3d18ZfV1Y/rh99drJ5fPb3cLbrF4x8uNv/5Yv3HTxdL0y0ef7LZXdVf3zz8v+tHO+oXPPbN893m8qIV/vfm4myxdN3i8fSl5X+dautVDTjnrfO+Wzz+cnOxvnt5frldLKUVv9v9fL7eF3/cnO2e1uLnm/Pz9jDAUhwfNkZfXz/efbXaPtlcLJamp+bbzZOnR1WfXO52l88OH/v+8vlhxf3HG4Cw3eLxg/3PO/UnH75/tVvtaKUUnjko3aklXv58u3q2PoabuqOOT88ddmh6dKweP8nr85Hhk59cbs/W2+82f2ujN6v86vJsXef1vtT/Pmj/vW/r2D5o/72/u7zz8Ore5uX6/Kf2xu7yVHGxjCXS9YMXfNIh279yWLFY3q4VPxy8xWz+sH9nXhyb+dPBC/Ujf9q/clgxNTN2ZXd5v/VGu3K9yEwaI97FlCWaIiWAnvv3H/xUx6m9TxEwx+JieZsPWJEcvTcxpRKj4sn55V+//vRuG/R54f4Pz/UGTz2Y/f5xqv1i+nXn4dXBh+48vDr41p2HV/vX7jy82r/5w8tnuhoY4Z/3P18+W72sHfrh5+nn90/Xu9Vi6bRjTzft152r5+tHu29Xu81l7cXXL549XG/r7+83j356uf/5c/355eWTVvnl5ZN93d/q3Xurs3urzQWrp1s8vru9vLp6utq0D07Fe5eNEs3RGuSq5T1Gf3V5tnm8WZ8tlo9X51frbvH437abs5eHxZ/3xTsPr+5eXm5nz392ttmtHrLed9sXfODzzcv12UG/x0/f226ebXabv6yvrpG8LzdXUNSR2rbiartdLP/8793i8vmOH3/vFo8/e7l+dLVYXrw4P+8Wj7+u5PmRQC2/3+yAo5W+fvHs3up8vduNNJIB+3r9cne99tM/fnfvyzsPFsvFv4w/u8XjTy9fPDxff/Li8eNxwr5d71abC0ar9fX+1eZv6x+uxvsPDot699v16nyxtDSut2tZfG9zSKDvj5uLs8u/fn/5nHU0Lz+Ylxu92j/wxRrq3NDgr+NKv/t00Zbz3dVud22Y7+x2dSOjf/c/We/+ul5fNEp9UNKh/Hx7+ez7y+eLpfTgzv2z1Q5CqIUHY4Gd7E4tyN+7xU9fXf5l/c3z1X++mDDip2/XjMth5eMvNk+efkkX2qak2LnaPXo6Du5P3z29/Otnf1lf7L7brXYvriYk/OnOi90laDA9+dX64sUnq20tgyR3HoFm0xuPv12vzr65OP95fOPxj5vd08sXuzlGjlj6xeqq4dhYM3/qz0d79VvjCKCYr+QIflw/1FW/uXjyKrYAzLh7vrq6akuC5yofMq94DqouzNKG0LV/gyxNZzozWK3ll1savWsGr/dsCEOY3jFDbE9Sn9qz/M7LXDor0mU/FH5n31mJg9Cg79q/QWQpxXaSbSfJDmKXKXb6/0HcUozt2r9B/FJs6tq/QcJSguvav0HiUlLu2r9B0tIa37V/g+SltaVr/wYpS+tt1/4N1ujDUnjfDFYOi3YpxXVSpJOYB+uWkmMnmaIbrK/txtBJMIMNS8mOO52EMFigakB6M9i0lBA68bYTZwabl+KCdjDHwZalpNBZE4FkcGYpznbadHSDk6W4BASdRDM4uxRnOgm548PO1XZyru8yVLa2E9zgwlIkdfq+D4OLSzEeCBRkl5ZiSqcvODu4vBTrdVytKYMr+1EPYfBmKcl02t0ig5clU6wwuzB4u7RWOh0figxVrkA6M3i/1CFm1qwffFjqmAKGlMHHJZ8pvitx8Kk2w0D5NPi81BZFwRx8WYpPnXZC7BCMAqEoNgRZio+d4uwQrI5DK7hlyeNvv0x+/B2WQcbfcWnT+Dst3fgzg9jtM2VpOg8ChyGapelirr9ZOkXqb7ukU7qshsgKElfqHRYR2KGvBwo51UKkDamNRF1JrrWiresilSHSvBZsGhLta8H4IQEABclmSECghRCHNK5h5jcBgd4xeUhAQCGHISkAIXRRhjSuZJeGNDYvbkhlWSHszJDNsgLIb1kG7Qa/7TLrMPCbJVypixkyWKnDQKEuX9ozQ2aZxLHA0q1QmSGzcHUYKLBsdbQ7M5RGuFymIJUiiadQiZcUyFdxekdipOBrwVsKlYqBVGYolYwV/XLSp5JQn/W3T/wu+tu6zgxizFJEOh+UKIgRxVBvO28HMbYShcL0yiDGLXPqRLGG2177ISF2jmJYWps7l7rgBzFxKcF0iQVaBjGQjKjPSuRTLE/bOddZS5Eh6SKjKoMobQ2sITuIQMakc9JZSlZHFWIUyyDQVn7zz4ZBIK7MJgPPdyVUiCs+C9TVSFdKJ8YMImkJgaM/kgeBuoJ3sXQCRS9LyYWGZRBrluKhWb6TCH2HjtVxyxTtUl8zqctlEOuW1mk7ZhBIaw5sCx3LXKCt4jolZZaHoy51K7nCDHWdbQHWyLIIS7LIYJkSSGL9N1gDbXJd+zdY5kToj/4brAmV3kLWB8ukZKckwAzWpCVrM0mX+HJeWlcqWJIGa8qStyDWYbBMCVjAy4MVWSbp6suDFat0ub09WCaF9nh5sOJ5s748WAm8WV8eLPMBRZfE24OVVKEdG81T57TVovCOzbLdQY4q0IPVKQHa2i5TUvsq2Q1Wp0Q7a50frPVTb7VlJqX1FqAtVHw/qjaN/a0v56nDvOzGNQxtM35w4zIGmdi3WMltF6/lupj1vrGDYz2P97U84wm0XFc1zwcZHAu7Pa7FPHEQWoS2VXYjyODNnjEygzeVJ+B1SnZ6k9KeT6Hk2U7YUQZvGI1Kveif18WstF/741nMlNmNgxs8qGPs1D8vppa5b+zgYZYAkR1Uy5a2Wve8OC0Zr/CLwsGj9EZ0k+NlLY2k3gxelNTXfgmknleoZ6PhBTN4yz5Tu2TblFSecfCwQ61LbGmeRdu6pF0GP1qX6n0oidUu1/twH61LDIHNY5fq42Xskz7tDEUa0IedaLHO2uDduOnpq053vTrfg3e68U4v6s47vcdw6KQNAZRsOMJXQpFpqmt5P/VAFMp+8ut9X/tT+zgEULK2oyCHEgG5ztoQShpL+i3te+3gEIp2nZKxQzRMge6KQzSVS+edaOizTtoQDT3WSRuiqdw6HANYWFeYghghX3WFtXLWLrdpGyJ7CvfrtA1RSZhiYG0SLKwYqF2KbC26wurnxOnn6CIAsrMoEraH63iM3ZL58SE2rh2qx5t1dbTOsauAjXU4rCICszZEUHLsH+9BsVr/wJIIyaorrJXbxt+mPyrNqv3ToVYuXVdYK9ftrSJItGXsHm3Dpdfeaamuz/ak2x+ioquDAkcS4c5ZY/q7Dkf9rajBzzQSwgr2kBohbMM+JLCugRhkSCBdK/JymPqrN/fTT/dSY28UhCE1BoeW9OHau3ozjzRQqUOGBrZxDTJkUymBTsqQGxFUIpHNNOfcaeBU2pMbNjZ4hzxio1KaDC423NM2psGGx6wUuTbRkJCHeZCttJEVnlRK2D65J4TcaeSodjezj44E2dghs49WtNeH81TSNiZaBJOrCEjzFNruUEdnyGyhI13nsxBJsK0S4wyNbEVebpxwA6lhYxuFITdsHFuqa5USIHFmbA3BbDfaVYECGVvfWBaZI2MbIpZFdnWQKt4MeU8ggdjpfgFM9dlGtdtGnd14PqgP72m2NqRHlPFVPz+jZK8Um+/od71S7Lr4huz1pNTwxI8Em376PThA5/fg1M9M4NAxP4FTHx7B4UNhhEYLbYTqFj1kjottXrQjYUIj/W4Aokadcz01tgEdcqgQtQEfcgAivYsoJChE+yIQNfqY9QzZ9oMc64ZWp23I9RxZJ23IepBshDPXg2SFdsixbu8VFYYcZ5ss7etxcv/VOj5jr/VEOb1aj5TjvdREISNAqe1+I0ipIvY4KKkidqXYOVXEHjuaKjM0jULiZDOesHOqiD0OBAfNehOMyhWx6/aXc5228dFcZ61uGjnrnNVdImflgRo4nDgVt+rkcuZsNFPHKNelP7ZRl/7Uhi791kRpe4++1ij1OAKNULf2G3Og21RuRLoVKlEcZ7PR6Gk+G40ePzNtQIrQjUbrh0oj0a2wp9CMW5mRaDC4QKTbXmzsUBqZrq2URqdbT8pIp9tUl0an62CWxjOMpcoxtJkujVC3m026pjgFUBw3GPj6bygQS/C4/huK1Rmse2ixfhkcx2pxvKpo7jkap6HYuAxGT+CeQ39a+sJBXTwfzUufOw6ePg7FlqWPXTCd+DwUZ5Ye8tlJkKE4WXrfKTB+KM4uveuUPKWhOLf0tuP4HZEk+KU3HRK3aIfiwtKVLhRkYUNxcelyF6uIsLi0dKnjiJxoIS9d6DgAJ1ooiE1i7CTFoXizdK6LqZNUhuJl6Synec7Gxdul0yOf5DAU75a2dJyFcxqK90ubVUhQEGsExFYpIC0cCpK02KXUSaGFhIQglU4KLXBY7TKYLEPxZWltlytfV1SO1unRLw0F4li6HJBJDQXSmLscOytuKAHpQZdzZyUMJfilxK4gvcpDCXogRBBmZSjQxdAV20oseRVPWL4Jf4fIpZb02FxiLUX4qq6kVgJlusI5H1ENMgNkEVbKUKKjJAgutKjAKNZqMdS7rt2NtVgFFyWmWqSTeSgx0y0x9JIiAtkO0YQWE6LPDlGplTSUJFpEHKZFRD4dIqFa1BESsa2IlLhDeoEIvEAZcycqa6CIdKHKC7WYapF2ucthvhOkHIw2lJEiYIShZKNFJAValFoEDO4iu1ahRy26WgQMPxQEctyl+xQRyHZi6T7FWIuAQTHVImC4oUAgeRgwKBYtOsBwQykVKmRPYofCISojfWzFOlbIm0SGApVEOg0YFOsMOsCgWGfQAQYyOEUocYBBUTFKpataVJQSD36XoXCk950gxTLI0hDaUQbhM2WEURCMsYxUvROPOELvu6UggY6dNYnnkWAikR7LoZbzWEYs04kvSMR5Pmk5AI6W2VUgO501gfsqTpUAPJTFLEvbirUoWlQ1BnftEmF6lZWIEadFYNG7XouA4rgbloh0gESLkWIEECSOCPCQJ47FrEXA0LtliTSPUUHEaI2+q9JBikJDkLPCmFoLkKg/CkOmrK6KW6UwZPC6tpXpE9wuMkioE2CrEK+VAVTZ3armKMCCyAjIoV8KTEV/pZE0zwGsoLbpJNO+g2NpZdrnEIZqSbUSlJEwVpliBh6HLLaVgcdxOKxyUeSTxqlGQOWkVZKbILyCeDYzGi5DllHcQMHFQN5pnzLwQOBpv9J7MZB42qdM+xB52mcHoH3v2CyqHJPx8J6tpJaBxwc2GoGeJ9r3kW1IMmNO+z6xSdWytp/Zwqq2ivHwhQ2ulhmPYNj+VP0WgScIm6OKRSPtB8tGWsu033hi1XExH8EvA++zodM+fHHqJHNwoP0Ql+yRWtb20jICD/e1vbwEd7RM/0NZouJC7qui8miQINYy7UVZIrPmvqe9qArDWqa96JaMDRo6T3vRL1Etapn2YljSlpZpL0YVdGuZ/sW0VKk9fIS2l5co/yirOD6WZUH7R5n2kkGrVMuML9uCAWDoHQ2yMXDA0QpaVJ6ZIVKKSIVfso3UCtpke2DR5KJqQjFsEEhZqbDaKoSFZqnQZlEL0CwV2ix6O9aoUWWbGDYKi2qRCmaWrQLU1QqaZbNQjSNPgFtsFyAjmk/k/IYNA2zkCaFZtgzQb1+BTB6tJU8w/mwbIJhWABgbh29woC0wbB2gmD4BYGwecGdawQCNMjitAFIVwrXxECBVeUgbD6X0uokwpmgsAF23kTamSsvZSBTzeAJI2UogS1oBpGwmlc6gaqQChqDNLRWqBFJKkDuBzKoaiKWtFWhdRk2QVqiGBr3D+A1UKewpqojIHfQUbZBiM2V0K6p6oGu5g6CqQoiByAlFqahCSBXoqYPAohGSAr6lDgKLTkgK3UodBBa1kKDHoaxaIsQ4ADSWq6RKywCsagi6BItHmZML8ETV2KAaQmvJ8xBgdEOohLUMfCiHFJkzSkZBO4SShf5AcFEQWcXt3GXgQ5YCxUedS/tI9xqqQ3BVRQSFL6aD4KIlUp0zZdpHugfFL9Ilxgt5ii4E6SC4aIksFJ6yto+GEZy2HQRXkO9B4SkzPghVdJm4DoIrbDCotorrqiLNLq2uGt9BcEU3GBDedxBcYYPRRQTTTzksrQP/ORJQjkt0mFqmfTYYXVKxg+AKG4ynvYhiVsRV2wXmFIIr3iytLrDUQXAFlTwUnDLto5PX9ZY7CLBwjtDlljsIsHCS8Ho66ALtc5bQ1Vc6CLBwmlAGRY+YIpwnVM5gusB4cKLQsu0gyMKZQsuu87TPqULLoYMgS+CMCmsROwiyBKSkyoB1EGQJSG/hRErngSegs4VxFFWtSgA+yq6DQAsqey37DgItAfi4HzsItATgo5xQZosE4KNcOgc8Efg4rZgOei0R+CjbDnItEfgo+84BTwQ+yqGDWEsEPsqpc4xHrOMHDqGklAh8prOwxsATgY+ydFBuicBH2XVW4QE+yqjNUbUCH+XYQbYFkYyWOQ1RBj7uw4lTZn5LB85AxCXV+cXWRYAnhXafNUg5tjIcCWXOi7wP10UZ/KOMchdDok9ePMaW7F/UIG/x+fnlaufsolucq/VYCN3iL4vln4v1XbEc/WJXLLQld8WWrjjTFVS8znbFua443xUXuuJiVxw0KnfFla540xUvXfG2K951xfuueA6P2IhA23JXfOlKMF0J0pVguxJch7lpQQceYlcCNDF3JZSuRNOVKF2JtivRdSViaBK6EmNXIrQ0dyWWriTTlSRdSbYryXUl+a4kjqKxKyl1JUGHS1ey6UqWrmTbley6kn1XMofU2BVoqdLv0pViuoIOutiuFNcVDFxK6ErhCJu6UqD7qldHDgMlNioJMWggVDJq0MVAgg1sg4H4Ggx2DGTXQGsNBNbAVxtIq4GeGj0VQ0kN5NPAJRsIp4FPMJBIAydsII4GimhgDwy00EAADVTPwMwa6J2ByBl4VwN5M2q4BCEzcAQGEmagWwb200CxDChn4DYNyGagSgY0M9V4hzfAKoNgxkBzTOANlcrBVBr2eIP4xcBTGnZ4Azdp2NoN0hbDnm6QtBhOJIYd3SBcMWzlBlbSsIcbDgiGzdsgQjHwkIat28CpG/ZsAw9uMCAwcJAGltqwVRsYSMMebWAdDbyngWE0cIFGt3ZlGtiWDcySYUM28I2GrdgU3oBtNOzEBobRMOdILWAa+IORA3NerZqYcwQR7Pf84cTPnCujiJRBlEOs9hDMufKHwpxXwwDmXO2lEA2wq/JHDSlogznn3M/2yB/aYM6VIawWWcy5soPCnHNYZ+/iD28w5xzC2Y8QZNAGc67MHwdqNhb+8AZzrgZbyudxPGZb4A9vMOcqDBTmXO3j1OJFmHM1cOMYC6HmD28w52php6fQahnGnCv7phZyyrfpoVAZtmpCwpwru6anM+XT9NwlKf/73//+9+5dWXOGX7LmrP4Iv+A/0iw2MdPeXvDgaNlca25NHzhyhKA4c3Y49nPYe4AIfhfXXUC+Wm1/Wm9nLiW1YvbJVjG5SXy/frm7c/EEs2sMUCnWm6Y3jIHeP988uVggFanl2fe5/fkl1uxRLYxXLzfXbcPv7HZ3qMd8+2zzl83V5vLiarEMQovcmX3wy9XD9ejmQntari14WtDyN48fX63V/wQ62yonsJ3CvXn005friyd4z5jeYOSsczC+qn3B3Pz4td35aG4+PTI2j2Hyg4+mhzqdv6GH/+ej6eE0Qb9yDj9ZbWfeQ5+stiNSqAk2/lgs0/N739U18el29dfqdlHL3zzf7V08aqF5edRCc/T45vnu02pvXz3IMIJnEeky+ub5rpEEOvHN893n6nHVHv18s3/ymhE9z2jl2WaH59hY/v7y8lxt6KmoXip3Ly92ly+2V8194c6uQXRENO/sdqxjpVO/QA7sr6QHLBdGoLoX0U9KOEX0Wvrs4uyz7fay+WyxuLWoj9PU5y8uHjXKwE2KM0JGsc0id/F6aQ/Tf4rtYZY+xdmUf7l+sr44m7vWAF2tnRFZPrSvHNsene34xEgZpgcZv4aN3eLxF7hDrK+OyHir/e756hHuANr25AE368Pk/tbqgHF67hCa6dGxenz0qGl97rjfs8q9X9EXmytwcg4PVXyvgRMN3R6fGxuuozM+Ota2B4+g4amvNhebZy+e/Z/19nLv7cGNA6dEpfLV8eXedv14vf23L/dP1/rZwNWKeTeBdF6772et/XT9+IvFMhgmbar5cbFMhzX3F6gCZo88aBX3VnPku7c6wC0an6r2LWvVdWfMe6uzg64zdvdWZ9edO++tzk74d95bnYHs9/fD02oeHNSwPTZ3JRrcPPqpOSvdWz2vXpT3G9GYKh4s0FAuHn/3aLteX3y+eqTUB/CgbLPhp8g6mKEtVfP5GN+arR8eobhfPdSMCFTRqj6zfQYhXcTeqy8NlUozSm6wKInhfFpB+3F0eeLJLyiI9mO33Tz/dP1o82x1fjX5FSlVbgyPnbiGWe/0gaPuad28fyCJVs46OJWn3aryKfAg1el3emsigfpO7U6ZnqY/4/fpTiU+l3+8uFhvv6V7PMlS089eLZZ/xpfo1i3+iL0lnPtv+Vt5LMdb/pZEStbwgG+P7EtaNa+35vDx6V427V42t7y55W95fVLbfkd//h3isl6drbds2Oo2pcM2lT7f7D4fkSY0pFFnLGZxuqO4pKP2aHWuLzP7///l5oLKkRe4u3o+L36/eTZxlCnnItkrxfjjs9WTNR+aCPzd1cXZ+frHp5urn9bbb1cXT5qbcq3/5PJlq6uzV2sVkpmn5p82l+ebi7G2eTDWR+9uto/Oj6l9u4WLKUDPNsD7cN2fvXx+f875jJUP5pUPTj05Vh48yYNfrV5+unmibvQg4Tfb3dPLu6tn6+2qUZ/rJ7Z7q7O35oLHjL3SBe/e6uxVRzYOU4zQjGJRHIlPuzsbv6Od8fpB7Dq5PkGrb3zxCSDAlVL02Yao0QRksbxtg+uRO5WQI+oIvFftYumkN3pZ6wJigYW67lubUl+85cmcChN26Ml/2/TOFOckxSg2e+j/zLP/dupt5HZEcoWobH97sWw3nXXBZadHyUO3f2tyH60rJSfnkA9qXIAxbIA1ri+Gu9kVf3C3+s3nEsS6ZJLLFvmORgMYfe9/WixjXwKDlL33DhXML4QLsJJ7kxkih5yJfh7GDrBOeqsPqKROb0+Q2tiHyMsI7ve3mI+jWyN4LZZAb7wNwZYUfRIfdWdrwQE09MCJaR5v637sUoheSvYpx3I99IB1zvQxWRusz5lV3t4mEMEpXBlvL5a3RXpJEr2xJnnkW4vH+2ACvTEVAUcMZHD3sQX6cHCbhse4BaZ31tUrhJgD4zXGMaC7FU0bnuIUXsMg9E7aWykbF2muRiswvddZMcZ4n6oM4yDuwUlQDyIhnIB2ioxwquEpUsKptlvkBIbreuQEqzhYQydYY32fXfISjY2JFfwrIykwdDeRFKZANb86ksLRaf71359iL/zHUfSFWn7d+Atvz1+/sTTvxV0frHslr6BH4ldxCyD/2+UWCEvDSpuHpdE1RlyaEwtsftSgH7XcjpbXWaw35K8+YyBq1KNfHLQv5NNfGrKGbS8uHl+en20untzaXDx/sbu1XT/iLWXZG/s+G9wqZHhDZuzXS8W/frQ+51CoB95XCLWbRPtVIaRad1+qNHXWwee71uOPWhSOiPHhBim+1R1VY+vUX4TWUeS9Xx8g/Mvjzzfbqxbr5MvV+It4W3Y8zD5bf7q5en6+msWq4XQ1nc8YQxX+72PHfHV59uXqYSv/gnD+9Wbq5+OZosUPZpKmE/8kPWhqjl/QV+wniXVb54jv3J9iHr3/KXqFduH1puhvH/QUMchKxj7yKTpSiMy0JUaFG59d7LYb5Epik4Mn/v7qxbO/LpY5CxRAS/aw+HKxzI6jhws2+ig+K5/Iey8JvWaiK8H3ko3lX7d4/NXqJXJhZeD1MNIExZVcfH25fTZKIMHopuiosaYef/fiGQAcmrtE7GCquUus1i7i1QgUywCLgUfA/8J3IRMqBMOj2DmMnAjUoFb/mPRipKVRETSSgmr2rcVzQk3rCIKTDa4oMWC9Ixrx6iB8IHoLogleV5P+EuN2xOPdW/1ljTTzlTEG375EAx7F+t7rwSaaVGI9fstiKd74PjlTctDDLCSGw7ppT6dMXDkwRfkZEdd7600KXrJG0Lov7XSU6hmIGHSOYEOwRO1WDEGbztHbmPQoZtut0t4qJjmfOfM+mG65+pZ12WP9OPJKdYkC6R8vVNkGN/h0dXb515kM6O7lFpHq6mxD8CrwbJrJh1vahsDW+GXEPNTJQKZ4Wp/O0uBmU6aD5Mdi5EmZvld1jyzR5cWTNfumKv730Rqr+ntz8Ssjsh2h05er3folvfkHzAyPjIzaq/i4V/a/N9E7TtuHAzGOy1zEPo0DN5uaouknZrwhc/ZqzpDx/ZJwdZVO1AkHbb7Zbp5sLiYabX3xSBWyCcHFeM2i5LeMTaOQt/71lpJIlXy/Nmb87gMSGAsJIQVvcnkrA/LVenVx69a/3rpVeu+dotEr8eQDQ493MBrf7c5ufbr+yy0dkdBHVb58tPgxRXfUiI6Hf8btb27Q8M3zQ1OIU+qZk4q+e6vtGkuj/1r8y3b9GFL498PGtbUa7rEt72r5quEf9RjyyeaCT3+23X5DQEruU/7mL+vt4/PLv3KcwSRpu+UY8mHwD2/9OP/dna8+u/XFH7/7nmX41o71X9299Xy7Ptsob/MPD/UM/KtJN4T9H2pYbg71b92+bX9evDnUvxcjRN3rfqWJ3n6SWCY3h/oxdjp2ozNri9/TTvRDnKLXP9RXzRVcWD3WJ8mut1b1ZDmoEmw84lsJvS+xKrucamnrmX6xDFmi9Kn4kEvE3rx9j+N+DL5Y3xsXpXib0Xi+h+N+tLZP0jRzeIXE2EczKupSSp0rvTWlnnCNw9Pd98HkUS+HL0DsU25jUfA2Eds73DTqlWyX+5JzHaxSsitd6k2W8YFiUhf6INa3N/Bq8X0IMr6Ca4nrs3ENDINfiEPZmOrEGHxJbJ+DG0HH9cQC6P4bqbO9jW5qheCSvQk2j83i8NDnFFSxrBcVScr+CddJH4pP0yuhk947u6+IXRXoTE8krXD7j2qFhP0ruZPeFj+Nh1ZIyfsnSgfeZGdrd2MpIxO5T8Hw22QoE6P6nvlEa1JIvWsrCJd8a00/TRaeM87EvoxjYEKXZVT6l1I8Xn2h9H5SRjvc6EyPsYFeGRFVKaBzu6hwsTcTFqn3p+tTFYdwEZ+UKa/jbNRjyYLfUwWeHaZPezzDU6sXNz1gcf/s4/4TBC0Jvd83YqWLuff7VyR30fahzS7oHTtGZ/YEgV16UTsEvXAaTH2aVXQ+9NPKNfhV+dDL/pvqhtjn/Ru1QvajYQqrvcyeIAYDg9zGNPv4DvxL3g3v/fDJVm5dYOJ2vvnbWq2cRinNiVszCYp37FU3jLhKkj5IBuIfMuKvPKQ2l07fhS52qcsd7nEdTmjEZepwPyNkR6du2x2RaQgD01lLrGSN3twRVCZ1Niu9vNHgvaHH0Q2zv1jeaPDeix78tZn9MOPzp5+jSl61c009P2rqeOjVvPrt5t30xuo5UpZUJZ3pJXemd96nGV/Zmd7bgDOx6X2UIrHDQNBhvdku15k+m2L9tOn7zvQFb+PpmUBNTsFODAtMr3FG/UwP2eEC96v6vz7FlEY+QetRTcnE9Xa2F2w5J+6sE+wfQ95XWLiVVM06leW3vXEuT+w4HL61ac+/840cg0sHrRhn4/yjITqZKuC9C/0dmUI+GqxPE6ScErw1Po5Mf/qYWe09ztQj1cgKc3WhF/742MUeG/3pki71ccY5Go5wM+azVjTdY72ky3NuVN3/pd+fgmpkgB4n6n5+sMEfuoepdr3bH4CorgfT6eIJGee1fs72E/bUCtfvkaN9Y0K48Qkznihqhe9nB7Ha7DTz4xN+Pmyd+H52MOUk9pFIox8+2d7aXe5W569mi0/en/HG8YY33icmvOGNb3jjd3lIuuGNb3jj93TCfG3eWM2blPGd8cAfFG8sfXQu7WXJne0dbMjIisDwxZyTGVkFRMs+pj2vGhEtZ1vi+ETqfC+ppIl10AqbZWJNqQih5L00sgu9oWL8BoyWidlO7G3oQh/TTLRcutg7a8LULBXiiUMTYWPyrPnYm+D89PFEheQ4fTwiUM/OTwAGfaJYM3K1PCEhFxnHBXh8cmkSZfOKFBvy2AUPxLnESbLPN0r2ggJBL55w3sZp9N0/B+8sriccU6qYNF6dZA5T+6uziP73Za3I4wBzqfLkgIm1vo+Hr8S9U1N9JfVERMsg2HQRku2Iu3acwaYHCmmremM7xynqqHpCUy6e2x/cVNViDtjmzBNhEs3zhM09EvjxUmjGAlemWRlxbfzoTG1SPxpHxKkfLX2avwJj/e/XLK1uXGMPsllfC4vwT52mejzgZXHJBa9OpXjGplLwgc3FGRNiph5jWxaNHk8lRxs9OlWMbZMUXGZTiCYbQuwde8biu+ptzkZyNvGaY6wzfXAiJZXoktrPjm6zi2Xuk+Ugm7P1gdB8x+mw39Av1mYXghfvnHFhdHsdHU9/WhCBc+4ZiyrjlWm0fx+/2NvW9bm4KCmEnDEinzxXcYs9McUHXrHio4toQmMI6lg5dn6njq/JhdT7GE2MLmcHJox+r7jFnkCU8fZieduF3sYYnXM2SqjBu6b83L3zzUuVuHcaLeM1nWJjiPXKxRAf+l07xV4H9J27xGJEfsIltlXj/RRdsb0xMUacfHXV/EqXWHD5xiX2d3WJJUrKpEZ9JLV84xJ72kfg3bjEqieH1NOX7nFagQHR9QX2MbvE3vv+D9WbolOZ5Xq7/Z9z5OP3XDCp/lG/m9JeHVferkvsvd0fnqwvtM8fjL8l0q43MM1EoFVtM8HZ6nBZnYlvnGLfXhDPG7HhjdjwQxMbVse+ajqr+MnP1xYbntjZ3p7jK6AR5t0Hk1yJOOPmEEgo5G0KAebdR0fA7eRFekvw6Wi9hM7bmHtCSUdOkeTackKA4ITnnA0EpS4WMYgrRSxh4sWZ0GebnZcYHbHgi++zsdaYbH3M0uXMyTPF7DwCStdlxD8+55A8gZpsl6T0wdgSknhPbgkrJZheJEsxOb81w9Ajb8jfy7nWtdhWgQOdhvuD+fHe+d4Hn8TZktXfup33qzY6hhCDekgrZxST9DEk47Cutpwjbpxr62Y+OZV+BM61r+IIX+k0+aE61wZvbPbBe+9tsdcknG/qXHuLS3nlV44MfBw3q9v1WPzdvI2TtykFCckU/1aGY+Zamy0xPT8eR9K3PhYHjrW+N4Goph/peEzeCocetepgO6qa3olbbQi/j1vtyJ1ILjb1gZgg0ZNcxNsovXMppYxviZB4D3P/hK7HOs1GYk1fjBdbpERLEqTYFxuJRmgLWXAc26JP3mWfbUFX5KX3GaVliEK+Gcloj3KxOYWUSYvhU29SwhMg4x0kpvTRxZSTEU9mHDGmTwKrkoJLpCc1GW2q5OJSseLelbnUZxqu7a257u43G6JY/H8Pn2yPJRCvfuKfWi7xfPdkrdk7/rtIJV7lTAB7ycEhdmC9ehAQA6cjva1RlwHf2XjjJvA2EpPcyDRuZBofmkxDvZB+s0xj9AqYBBljxRu5CUzsgvPG+R47f3GBxFY5BzSYORljHRkYUyI+s1isT1LnYnA9cg/nXdKU3Cb7nsjDyBQ0NVuw+AZKhMPwpCPDYTP5YLMNRpPb4SpcXLHRJR+j7bIrvRpaOReit0FlF1l8LORvS9F30ePBmRz29iaL7WLKtk+GM30qyb21UGET46j0+N2HQ5km4oZvWzz+5uH/XT/a/fHTFrZrihemjp/vlG875tlutEazSG5qwH5c/vCiRbD332iN3nVIjxsO64bD+tA4LDbRD5fDctn5PngJUdOC5pRMn3LGbS4X6WLCtFdicd6LF9e5aDKsDcyTC4akxj72kVAwyVkhNavPuTculWKxEyfJZ/a9NcG4BP/mu2IIz5KDM/hnOttlk/ucUpRsQkZAEwshMIIj7qyllSgWBzZnxBlHfmDiy/Q4gMZig3HVkX3ak9khP45IIjcslsYKr5ZAv8hivfW4dp/e+f7Ord32xe7prSYMGkNr/NMzWP/txF8cnd7TxvCKNLKvFwb9JlL9+0gmcBOpXind7JTUMj58OHEHX9uXz0YMTCYOK1d3gtE0Z17USPUhxd5JtIVUWhjw6ZOErishOPJ5heS9c56heMfhMKa9P+Fi4YlHIJjkC5IqOCeXSmdD7AQuihTv1nUpd95//BKl39LljyVCwcOHt642Ty5W5/+Qr+C8+s9k7nvDV9xkwFGB3BSDv6Wy+eWIxjdimxuxzXvizl+br8hk5ZyJboIzZFpscYLydIccNwG/Tbw2iSlk8darXMVi6UNSV8ZkQiYtynTLYv4ZTSCRSo6p1Iyn74nhiD73uUQNlkSI267YPmdpuVYsNjTEpR1jJVlUb6SBmaJWRdfhzz15PucSiS4whXY1cDBIbKY4pkVDysY4c9r2BOTKkWheKOvqyHIRFSykWdxQApr2Jdi9m/pHzwBZeNAxhhRW1s6U2JsWvyzA7GmE1tG33LrOF+mnmASpkIUo9mEf8YrgUiThmAaSYLPzoK/WdYHwCNMDKs/r7T7Iq5p/k+9nuqST0ot0MfRm9mK6cWvfHbqw32R8PuHzrC5+JzIE/4Jb+4k80IcJn70nQuA7cmt3uc/q+I5nfFGnxN0+ibLzvVBPtJBr6Z4Tsb4JAVGIaJLYBX7YzZ3afTIfulM72xFe7TaUbFyg/6Pb+ps7tTvpo4ZjRE2QNLzB6LX+E6miNebGmK1Z84WPtzWZtjj8VSzhDpIY3FJnXu0xNed08mh7QtK/plf79XTO7zbV83VA36dXO177P7yseTyZ3DHRM6P5K/3YVeayefTTGEZ08+inn6sU/8vLJ63yy8sn+7q/1bv3Vmf3VhuS/cD1391eXl09XW22R8V7l00lMHf+hec6yIeLiObybPN4sz5bLB+vzq/W3eLxv203Zy8Pi7PEo3ceXt29vNzOnv/sbLNbPTxfL8YcxZ9vXq7P7lw9Xz/afbvabS7335pnLz7yLPqlDH+HZjufvVw/+nUZAg/f36d2rtm2SDKyWOLHruUbP/bfw4+dhXXgx65Lao68yBMPkPetS7M++eY+iqvXMR63vy6d87PVbrvRLIIjtl3TkP2ukiynOva367l+I8r6WERZH17a7TeTZ/3DxADvUZl5oyn759GUHUu0bMilb1IroiGOsilii6U+ljE7DNmT9xIt632UPktxxmVzmOzZxhis9Na6VKI1bIrvWKKllB//Fey6vWsJbLKNSDqcN020EchwFPqyz3kUCGXpigbV1CsSUtOkuI8QSAT6nGa5mKhwpJOdLmRXJFaaKvQZQeqND31M/MdTECtkjsaBu4lRcIM3SGP2YhWNY2+C5YYxhVj4xiAg0x828aP9z4Xe2JYd2diYOtcH2WfnEWAVcQTGx0p9+lJv+E5Z3/bTpwxx0qdAoiXnQjIsvl5HhrEyvQ1lFktdISTHDB/UMP7TfxqYUxOIAWcRHSOh0oufhjkRlz/iQzxeCnOaJbQivKPpjSQ6ZYyL9eOthfl/hNwB+wCTkWxaeTrioWBNvZgyyRENyQWkzDIl6QROTbUpoINjT+fNmb7EcQFVmafpk5tJ4EwmhCnZh9ql0s4S4n7OdaYMZnF00et/QJlp/Ez72ac0y/mkQlJvxlxaXMhTvfHjrKlDRB9lj53CdImkOpqtKccQaIttKrXhsfU+JJkyKCAeBBFmeKCI4BLK9PFynSVt976GRRKsnS0SXR22wsHgzpqseDQ17+OI41wMjqRD8aQhsdt+OHX+EM6M4DBStjcyRZWtc5pkQoIRv8ZlkipWj8tkgsXl+VLVpU8EOIbU0iNrZsmwEJyqneR08YKLcSZbpcZYzXRRl3lbRHVl9i6oGWa7tL3akEl+tqC0PyKaJ8PH/bBXhLBhCoVLuFJDar5KVxpuh4px4yyA587MMoXVpV4ws9hffCckxZxURhG1XgyHswBbykw2rVhu9wn3NLUCvVewD2jSSElsmucB0zlJU9DesS3jReF3WUmcLi9b5vPAiynYvVRcVw62qLQ/dp9VoH23Zj7JSnBCpUjO6Xy1/3jto8mz0a0zEeM89RhQIemnLbF1qdW2xm5KkRmu1gbrBmCL9sl57bv12nztr+Q5HvFSkLkYXxeYs3VqmeGGxWQvqR8N2lJbg42ySt22bN01guMRH5paSC9dgJJ16EYCMqf+kudzpk2MLdU9SCoOtnfrqxUHxOvXnZeZOkK7LEH/syeKvdg63W1MRHfSqReVsLXRrk3gK2Z65qqtJt0LDPkb9+SyUQZTl4dSAXF1i3OiI9oGqP3HZO3gYckqKZgCVLdLqcMsmV4dmjZyirQm5+N1NCKAiE5Z5SMmQqFchUEQCipVUJSWn0IoncrxqYMeVyowLsNGE9qecACf7uv75uuT0emCD6kOTx35Mimt9AK+9v02UQ3aow7VYRgrK1mqqCjR6wg0fiNPeXD2328cTgO8UfA6RON/4EsqpZnW3shjRctMCpF5DIEB5pGmFf4Rmdpoa93Y2LhPzD7emqo9qigvIaLETL4li+QCrV1LsVovakZycYLk+gpgI3at021sG7VolRAN6b2fhemuTtxvwUT+SPL5YcVUstJnZVichovgUPMLMZVsgfWCGMRQctCIy/L1p3dRR7Q9UKIpzgaNdzveiqXyRCn5HPWYc9+2txq7IsUk5zU+7IPplqtvWZe96K1j+eAfLza7GtLou6ers8u/ztza7l5uL9bbb1dnmxejsfzzljb14Za2yUy8eoYkmFAPV5RfN3JO9VBZvxxD52iAxNeJqYTaYK8FNVg3YM3266TbR+j0W+IG1RCbm4snc0nlK8MGwQKpwuMDih2Uiws2JjHe+/y2gynVg//HFUwpWKtxqnMOMb29aEovb/3rrVsYSYRfs0BUnrdfH2Ox2k9Tul/VVw/qf2aRU1hKrzbwRKH15ebZZvf56tHucluN5lCsfbPdPNlcTEZ072o0fmY0bpU++f/mozFGl1L0cH0yGobn1SRkJK4fSOSxt79YxgFRDHGo5j9iDJn0pyfibWnQrXcfeaukdxR567tHq/P1YyUfLRgiQUl/XixBitAXK8nnFNMkwv1ZM0okIe9gSsFHtVSp9oq8lkga0GfyTfqsqog72+3q58UrIwhNIlhJkWTeenwTkzrJ3nMc1iuH2AXn+n2GypA6SWWegzx2scyt0chcXOZJycnhaZsNTbu62M+PwkbIgzRL6qJJPJFo9qHz/V7KydVZybkP43mI0KaaG2h/KccOlw5jziEBLr9UUYH+FYlpn/pGcieS+pnEzHeuT/NcjRwdVACjR456klABb/uoDa5l+9CrdMH5McF8KSnELvfVcLNdyJnquUulC/qnQdi3L2fp5+de17mS+0n8lXMXfJ/3p25BgqfShsCRSaWG9STeJNHTf3w6SPfeEXh2L9+yNVH9TBSXu9iHmbSjSi5opMpWVAjT1z5MjTg7md1yYfbY+5mYC6cZj5B7vCKZfvZST16RPnJO9vypR/lpzOmaZdrH19UNxx0kEyqdj33e94xM22aOd0xz2Uuw66m6tobwpjapU6OiDBoXjh5d7meCKRUa2v4w+VIpvZuy2rIAUj+TlSs+zUTQo1ztoK3ZiILgkzCat0OPl/feDKp+wfqDGiLakNVrvAgO1qsgqaKyypn0ULzvHvnJ9gupdLiRz8qFhXG4TuMcCRmLkFogXL2IS2YR9u4vVdCMQHFxjK8Lq8nzVNkwP85fb9b2BBEkI/DsYqDm41Q6T2KqfVtk7VURoGJulZSMg96J6S351w5FLdgwW/K2+c4fJq8q2FvP5WxkLZuLEbVjVaoxDjOUKHZF0bd3yP6dwqQ5irtymNMrg+IzkTyLBhBbCmaQtDOQ0hmdUMqi2YnJLF//eKy6U+ysnaMDU1UJhooq69eALvYz7Vn9YAJX5mK4Wq1S3Ew7DiIxl2UCR+rNXP6shuMNeMbBdaFP1Kqk7PpyUBKj0qDEGBHxURPb7S/W7+F8N2pRJ5WumP36bzBDuHQBREZSlzs/kgXZD4m+yFHWNSS2E2Ugq/E1fZeK9FQMVsFXAgY1QY7kbFf224JezAzLT7cqRbE53uuYId7SYdA/KqXa/0qhc6zB+YICRLpkmXcVgx9PnSGtW6NwVYpGIyqHs6SpJu3fdSRg5GzqAFWxui4ilXO35uo+Wz/p2JCESVbhqifX4UFGZCDQ7+hwab8UDAVoD5Wu0oRKjsCi7JqoB8dJ1OaVlOm7lWbratM/wKWAGNIk7tMbcpF9e5oe7Ts9q/3Rlynuv6UQRfIOWtsFoq9PdKURuTqj9apr6wQxIZlfr6oFRUIdgOvTE3Vd+M7nOYURuTHgvzHg3351ebauMov7JKCzvU+klYsxeBs8QjIs+CHKzppEQN3gAgIKZKq3EbVaHg9R8wnwOO5QfZBirHG+ZFfz1R1a8OcEAZNgJTqpqdPG1HM/LZY5+t5F7qZMID1aazb4i6XzuS/ocQIR87Pmwjswsrcp+T5JlJC9k4j8ZW7BzxGkFO4GIgTO7tKfSCqElHJOMVvNendowR/n9vuowd4kLZ2PqbfKdYrR7H1/mnka2NgH5TasAjneYo6Obu0Tx2GizxSWJC6TBDUUFcnOLPhPzfB4G5m3MTZHDD2ykeiKpibcN4CZvvRklrMSTLImCxMwWunr3WuYMt5dLG+rAU6IPjibQ0qqWD+w4kcfo1dK0WcbFAH3qetO3b7/w/MmeXdVYOhcCDEHndrpnh/5GmtdEMT8P45vSUuHRz43FQt/0e5MZwzvU3VjvfPwauZycAqaQ1v+Vzzx6qbvPLx6des/vHy2uagzrMEnYzSJFMDVjJ9bgkse2JslZ9F8Es3Y/7b05cRLZLSTPpjgXSEmt8lRPRSfrnerxdIh5rr3dNN+HVjEc+vrF88errcLtX2+8QT4XTPaVfnt3hNAyzeeAO/TEwDSerzE2uZ5vCjb5nltRbbN83g5HmsJ360XAXLe13Ah+O75+eZirR4nP1xs/vPFmsCoUIVPNrurhWYNbOrI2ZMjhp5fPvnDo6eb//W/7P/8H/zerV5oirjXVpxA3T/fnI9xP2mX4ne7n/GF4e5Xq+1P6+1MgVorrj1RdSy88en6HKp3m981UuepKa353k7R06+fN0L5v6926+d7p6EvNle7yyfb1TNEqP+x+3lPKb6QzxntU+MXrFLYqtCtI6gC+LbEP8gB/PrR+vz8arEUzWL4dr0x6Pw/Tah+pva3YNhNgkH0pauHN8H4rxZLdFl/7xY3Hhr/RB4aTkmD6uLqZqo/XzvD4Bh8/y1H4zcWPJvMjSDFH3JsVUd65/M15yE2oufKD+AignW8JOKr8YPAI/ywptWgKa01jngn3CL5nv5QIRr5zhIyN2+caTVuvOW8ipY95tXtmWLr215CrUEgq6/7GNqtgsCv9yZgEKY/fPtySA3AaEhT0HsTXa4PxyRq3WcSOZS4lbwqubzRbEaElMEJXn+EaglocvH1M4XYIugfsoZV8WI4aAs/Etkeey9iCSbCj6RwiiWqDRW2QSXO1WETp/4i8L71+xLauyHV/kt0SaGUhCkmj6akugov5I/TmiKlPlNK0Hmxpj1jJXgdGvUq4mHrrKs1nhwP6vNhSNqE8wceFvyIxeqXbSpSazJB9FR3VufZifHaKWfbYDlXv+K8J0lU712IKIV772KugLosddRcqUiik8Aj3po6qZ6UmFR4TTXVex9tro8Q4lhvlTqcweQKXLBOddRqr6k9CSFWGEJyqJWJfoxGVXofsTvUHzbU+Y/eqTCYfBUVZUnmqU0hzdfGk7SRTa6ousinUCqSoKDWJrIxqkbx2WrIwd5nTyREfsTaVM6pDkkxJeksYB2pg1RCCApgSbkwskHtKv/975CPf9uunj+tuxas3HMkD0rt7l2e/7xYHlsBVpaXe688YigXdVoe8gCxhjcSfTCC10FESvKJSrhstjYHEoaB1tiHqrWnmFxMSNl6XzzGnhxnSKOaCndKcVbKNTu0XwUmobZJmTqJbRqYFu+n5B0ZyCYwHcsea2ExMRbqK5heCpG1M3OlUsN3AGcV6l6D0yEX9CI5aoCNOpw+lBRNCFlTtCGwq3AmjrEmeZK9aSwPhZNFlyUhPiyWfeXoRPTrxtPFQNSP4/HUyOSlOAnVyFanPZB21lpqHcbRI5wl2CjGGWuiVzvYBmcIAVFmANLr9oe/Dk5ryql5D8GzHWXjVYjX4MySnfigeWK1XsfTSohQj4isOIC1FU5cXbwvFqe2/GbDaXOJiP6PhzMFXS4p+MzRuIKZnIhI8C5am2VaRTZadtZorTdJR1/BDCHlFLxSeevdG8IZ8qllJEa8I3GgC4ImosJZRCz45gjKo/JkHU4n3knxeO8Yo9UKZzSE9jE2xRKdf9PxtFWVcjSeunO6IlbKHj1x7ArOOZuMJ91gw06XSEeIi5WLoaiQAmSO3tiUJYj1JfvwhsMpuCZcm3YJErOkQD7miSqJ9ZiFiStGqmZGR9MbCSaIJJ+itzMwURZplsXonbwZmPjvnVhEzG00NlljlEjqrItL5E22xkRrHPBUMJ2kwHSkVIKbaHxkjwg2FrJVmzdc65KNkxOjWXzySYKQN2pEToZYiqbQNkEDGlUwo7XwlzljcYA8vk06eaJ9tNEmF950NCPpq65NOkm4i8noozVJdR3NGIOLIQUjYt1E4X2xPkj0ln0rQlArmM6EYEuKNnhv3xRMX8wJ3CTJOBquYoKG6KpgZutyCnh9ZKueHDqaQdOSZ3GssD1FiibCcnEgMFbecAmJO7mx28TcmmJSnC2hUoQYVSmWEtV7voIZjLfRRkdQSFXC6GiGhFOS8xbi7K5nFf5V+5BIyidw00EkozESrAZr1dG0EpMmaQULPUhYwUzJivHWiMCOjpMeLIAa56y1Pr0hfRejOuVjuqlMuQ8Jr/RxBaHMyzEIfmI+TQs9GsmBZCTBZ4JZjptlLBYsTtjPufhm9AhV2QnUdOiNS/KSgi6VOpjBO5fFmYh74LRZgoHB+OxLMVbjBemcc+60DlbeOFPecDANyaivL3RvXbE5EIktMD4VzOTFZktOFsPBss159DaRF4ZVFyzKfwXT2eiMt8nEkuKbssYmhFNkk3XuacdV3qyCWbwPsHI2Bb/f0mPIRrzNJOD2GrFCwbS4Wjmb4QNcfMMtHT/SE2QzSCiRaSfK6ziaToIUSxBVToJ73IzKoNjoWXLjYEow2SW2yxzkTRk5Y/0pqhlizCFY9kU169bBdDBwuDLGjLn2NOds5FZcKE7sfg8yATa6EGEvGf+mjJwxJxm5CIVxhCkouqQrnCG5kmHRrXAeHnET/i2H4EMJ+MYrbt4m5oZNJWTH4QMu+03WuvSct04QzpgQrqSQc9aTWAU05ZCjg+EEgD2gpKzE/09cVov7T2tQQ2+Kz8laZVPfEMxmu3PEcOp53zsnmcN9W+uuFOw7kBNh67IH07gSgs2RM+l+PK3TkBVRDOPwRjRJSL95inSy0UggHEpUEyMdTm8N5CuwLHza42coDsLpA4ZKezidzdGYEp1zKdg3HU9/kuXMgOiMWBfVNqfCyXnHRMvWDoOnu2X0rmDs4l1JeVxFt7GYzDmk4J11xps32oikLy6TG/141gsk1YgtmDONs+6DSbFYl8SmvD+sR8cpXgxBV5yaPVXk9BYQc/bWBV/eiJcDzpOccfFCnvXMimV662hGE01CcmMcnuDjeFofCT4CESYow7jafXKG3GEc2oN7IxqPibYN4NzRKioJ87eEtdCexvtkbDYxZyTFeyFNFDLFRutiMkj8TsCZXHhTonTN2KDa8gg5aNlfagjaOp6ppCxweMUm5S4VP41JeaKaGH8R+NUJLKjLPuI39fhPq/NP1k+qNuJPq/PPLs6amn/95O7lWPjs4qwV3npQw+eoG17DHqHKF1+hTRfjIHFNnT49OmrTp4q3ZX9gVLJ5YIKAIHm0SmhGCGwOzSph87dmtvD188vNBQYUYM39xfLPty179IxeMikPmjZjLunSBFcvLh7hFP7rHLAn9zjNKjwzV1B3ucV1ddK1mF573dBbn/+H5825/TVwgEgE+Gy+Cg3Y/d7SFKvwuZrMEi46FjWYRX8ymswiN+Rkha1jNWLFZNZiFOjFpURYqmph20xmM0xujvA3kLYpDgEyyBCUX49RifYUhwBOzJuIGAsOjrfGYAPJZ89CTuKd0zPQPg5Bi17Azo9qZoxL2tyNu8VvjUMwTlVbZzs1VhldZS+eYJrD15mifVvqZFv9aC2wUJ6Z7VCc5alS7X/LwPzl5cWT9dVusVRbmQ8hDsGrzHK0j68xAOPotP6PxffmZx44e3qLcDGZUA4dZ28cBm4cBo4cBqR3KcYixmd0p+oAIIul70tCoAqXUypbB/GTPtjo8SqOHo0RpNIuloRhjAZVZkacDRE7sOi/LSnG3mIwHvguVHx0CFCLd2NCb4XbWQgnuL+9WKJJ6YNqLaOq1tUdYDRpJ6a9+NhblG0uxKBB0ef+AsUEUmzEXHm5+rqatGPP7jJqWewCJGeVmB76C5BOoLqVQYyVA/hhb+TP8h6Li6WV3Dz0UK04mvrTwTi8E38B6a3Y6Cy8YXCRwRv9Aao5//H8jncXy9jnGIJ3yeUSXGzW/OPQ7tRbIPTZE00u2VKSitxHd4CfTqLJeFfPxoQG1PCDCHtCzcYws793vlnuiylkzTsI+X/K9H7vKlD5qRhjLkYU5d6tq8B1UN+DowDeAIIE2ibjyCKo4f7xEwjESnM+SJJQVOnZ3ATw3BHvrfcqymKFtgwBoXcBkUZx6kLPKv2VXgJwRjdeAr+rl0CNebP3EtDyjZfAaQtskP8tHRakSno5L5xYlOyY11dk2zCPl2PbMI/X4vt1EVidv4dD+Ze3Hr3Y/kWt3t/eNLDlHhzLOWb/Q8+A6VCuHlH3fyEUiun2loZz4yOskQ7MfKg4sKehYn6c14oDCxGeOLDFwJLuwOiBCuzyJusCKg70+Gp8N9eYU3Ggm8ZC7kALTMWBvpWKA80m9nEHOkQqMKEbtXWUD/Riai4310Cpldxc14O924FWhYoDBQaGcAeaAioOJPKCYLC4vewbM7YDITMVB9JcwVNwLjil4kBCScWBKBAzt7nMTYUoD/4xjoTeEYC9kCLPYjUReiewmTlmUEfjkxuCEsZMbIDAA4ZY0pCT7IvVCkOY3Ggszns4lOeCUSEGfV7LKflExMyQc/E8EL2JnPyzcUR+6W2IWY3MUEABhcXkxCWEhcXSqvVijSuSkTNQRp3mSvCSjdcKi0GKI8RziIBlxWeXi1Z69XM3wdInr+YgXcCcA7F4NOTT4gnJPrtAhMSQJFGRvLGC+V7BRTgQFzOzRAQ5pT6B0NUK/I3UJ1yRIjaiUfH6io0ZdRc2f9qIoBnzFlPaDKAYOsTsYzI5GSpMTsF5413AXJacdziGpiicXLQionPyqcSYglaEkmxKqgbSb5pgA70rMeeoT2DViVOqoPxT/HjLgrojK7pf8itSdu9VgtDTNDeq+AgiPPOeovjj5mz3tIoqf8EZq25+v4Hq3rgTna/nkjAYt5kobEzZrptadSc6wV1Uf7UTjP2NM9EH6Ex0gg9UJ+0Tx6/3P383uWBetQQROenJGv6yrTh+PdyginmP6XpeO7uxipaJFPhb/IxOYOPkcnQdg7++3D6r8U1rW3cvL3aXL7ZtXL578QwHJx2jvRrpY3IxUveGqf+TD9bq5ebZi2f4E4uaCew799Z1ZP/ddaTHZ6dRRXrMMN+oSH8PFam6T0XikhRbkpr2VJGHS9jMSCLkjNrwVJmHx1o8lRy8x4dpryTAoS4Fk6o6c9KPeiz4nU3Y2lk1RN7rR1twRVd8zjVl7KgfDQW7yGRTTBlPM5ppgV/+6fSj2O9/COrR2fHiY9GI+pjUaK8lAbjRiKI/MD3L8tvNkxuN6LFG9BU50E9UQ+xsbrF+RS221HjkdPWhSvRdpkD3JfelZFdSSUhoZirPn4gqnXubsisYk+aDu4vl7dLjVlFKTs6h8KsKzlEliMrvbaZAfycK0XebAj0ERgghFznx1GJnVHn+Yn71qhAVMc5HaJIr1cVlpg69nlj8tSOnTUEjp+ho71Ydeh3U96AOnek/9eeJLOjU/0qtpr5ykwX9ogYq+l2yoEc1MdtrNbV8o9V8n1pNVgE8Pf/VfU0r2Mn48X6Vk287C/rdy+32D7vL3epcQ5SNiAbSzfSSv7eMXF133m7IrXu7PzxZX2in/2kib02COsVUFa0qir5/Sd2D1cvNdTPtO1S+YukCfFPc/Kwvq0HNDqvaD3KuVMQIfM3Wt69BVLXqddQaH85c3Qi+ddJO6J6mOaqzpfFkqftwBd+VUv5m0fcozb0m730tKfd/Lf5ltd0ulovPzy9Xu0iU1hqoyvu3FGLrSBP8geVmJJ864hQbMUhQS9163K3x72MIxG1oIj9LLgJ92kqswcgnmd9NbsbZYjzpE/HB5GY85p8+prSMrmXxkJR8CNfDP/2WXJWfXey2m/XVrX+9Jd7funXr40rLGEuwpgRHhI8TYTN+y4B8tV5d1MR74vtA8L6PRjr8zoZD0+7dDMdiOSUi1MyMt2LvPq5Ule9gvRykZvzoR2RyPv39cjNKzu87OaOejTRR40Fw15+rhcQ/SLioLKwGNCWVliMwaMihuJCT5j+z0RLNSn3Ag6YtskmM1WBhCUcQ01uH84t3Sb2qyLvkjXPGlSTWORN9d5v8RSkTPSA7gz62u+16S1QL4kzgIZ7WtzXxW3Q+SbZGHbsseZBsxuDRGSLfEeHU9EiWjWRbSOmiX88m42gSjMazPNETTY+VkHfjDkUcQT7tBFNWH4vxLmuKLOsIvUnQKZOCJpJyxVsnUggf5aOrfbFE4EkRqHSYjJHirPMhEajNC82ZiJ7YOgxaiYwB4CkVGw2RGI3RLPPWmpDwQ8ukAbd8PeGnpBGaYjaW4bw+CcddEQVLMFVlmHBf0/xXoRCqI6bifGAiDEaqLuqQFmIj8Z6atBJ4xNOdqLBjXVtcDDbEUILm9TLWY74aovfiomP8iBBrnYsu5SgiIIOJkViJ4jAmJw4tc5/xqopiigNPGPhrSHQ8Fbx21J/aSXKsBB+JSUpgVVLcMznEFo0lWW2QjPUmG5uSBEcKMeK6RJOCEM/HEN0JCYLLPqRYRB1UAZ64nVgASAneKDoWjQFiYyImR9L5sVGIrIUlbyBwxKkFcB2Nrk0FHTzqS+0gUXEIk2G9J6CtAg+D5H3OOWlWNiMme0N6b0e4On3IEcjVJC+ZgIM6+4raYGQgHI+mDUtYVot1Pov1EtrsJC8m8G7UfGPHy5cgI4eL4DoenZyH475oB22xxeFNVryJtTvBZOLalpKx9G6TI04kB8lFkmgGNKzGSWheIqbhtEgQYBtdiESz1Y87KZr/jOBhxLM9RXkgRocL+NQioO4EGs1nwpzsT4XDiNqjk/mK1FCmJzpvtj47gr9ila4NFE9UQyvZWWzEmWy1pg9SvMHERXtuI4hKgGbvwfBX0M5TROjUQj5aDLR6AqGOJ4VYM0e9ql0NtniTJFnCFdVVE9EtErQvl5CULEgkSRYOA7nEpGuSFEtecoiRkJ6Ksdep/wkqeooQXV/LpxaEu4ZUp2bluEO1l07AzpBzFGeYKXHGFONSVov8ovPpbYHUBU+AIaWkJ3ava5vACTJ6nRJdX8zHC+I6QtWEg0eTwTwedabiMRMhwZeSsovaH7aZYCEohPkCi4+2XioON7ATOwDZCo/IqG4dh5To+nK+vhgUR65j0omZOO5MnUONDIRDSilWdCs/YhxObL+Qi6P96zYZWQ83gfrUAR1lnI+oEZTncCmz+x4tA6bxGI/4/NFc8Pmj/rQ0jFgmfvdou3lO6BfVLX2xefL0HAubu5cXF+tHu30+mc83L9dnVab+eHV+tf77/wPY1TtD').then(json => {\n",
       "   const obj = Core.parse(json);\n",
       "   Core.draw('root_plot_1779221492071', obj, '');\n",
       "});\n",
       "\n",
       "      }\n",
       "      const servers = ['/static/', 'https://root.cern/js/7.11.0/', 'https://jsroot.gsi.de/7.11.0/'],\n",
       "            path = 'build/jsroot';\n",
       "      if (typeof JSROOT !== 'undefined')\n",
       "         execCode(JSROOT);\n",
       "      else if (typeof requirejs !== 'undefined') {\n",
       "         servers.forEach((s,i) => { servers[i] = s + path; });\n",
       "         requirejs.config({ paths: { 'jsroot' : servers } })(['jsroot'],  execCode);\n",
       "      } else {\n",
       "         const config = document.getElementById('jupyter-config-data');\n",
       "         if (config)\n",
       "            servers[0] = (JSON.parse(config.innerHTML || '{}')?.baseUrl || '/') + 'static/';\n",
       "         else\n",
       "            servers.shift();\n",
       "         function loadJsroot() {\n",
       "            return !servers.length ? 0 : import(servers.shift() + path + '.js').catch(loadJsroot).then(() => execCode(JSROOT));\n",
       "         }\n",
       "         loadJsroot();\n",
       "      }\n",
       "   }\n",
       "   process_root_plot_1779221492071();\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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
}