{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "170b2646",
   "metadata": {},
   "source": [
    "# multidimSampling\n",
    "Example of random number generation by sampling a multi-dim distribution using the DistSampler class.\n",
    "\n",
    "NOTE: This tutorial must be run with ACLIC\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "**Author:** Lorenzo Moneta  \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:26 PM.</small></i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fddac29f",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:11.356653Z",
     "iopub.status.busy": "2026-05-19T20:26:11.356525Z",
     "iopub.status.idle": "2026-05-19T20:26:11.680504Z",
     "shell.execute_reply": "2026-05-19T20:26:11.679706Z"
    }
   },
   "outputs": [],
   "source": [
    "using namespace ROOT::Math;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b850901c",
   "metadata": {},
   "source": [
    " f\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8aa5b520",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:11.682473Z",
     "iopub.status.busy": "2026-05-19T20:26:11.682352Z",
     "iopub.status.idle": "2026-05-19T20:26:11.698236Z",
     "shell.execute_reply": "2026-05-19T20:26:11.697703Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_53:1:10: error: use of undeclared identifier 'a'\n",
      "unction (a 4d gaussian)\n",
      "         ^\n"
     ]
    }
   ],
   "source": [
    "%%cpp -d\n",
    "unction (a 4d gaussian)\n",
    "#include \"TMath.h\"\n",
    "#include \"TF2.h\"\n",
    "#include \"TStopwatch.h\"\n",
    "#include \"Math/DistSampler.h\"\n",
    "#include \"Math/DistSamplerOptions.h\"\n",
    "#include \"Math/MinimizerOptions.h\"\n",
    "#include \"Math/Factory.h\"\n",
    "\n",
    "#include \"TKDTreeBinning.h\"\n",
    "\n",
    "#include \"TTree.h\"\n",
    "#include \"TFile.h\"\n",
    "#include \"TMatrixDSym.h\"\n",
    "#include \"TVectorD.h\"\n",
    "#include \"TCanvas.h\"\n",
    "#include <cmath>\n",
    "\n",
    "// Gauss ND function\n",
    "// make a class in order to avoid constructing the\n",
    "// matrices for every call\n",
    "// This however requires that  the code must be compiled with ACLIC\n",
    "\n",
    "bool debug = false;\n",
    "\n",
    "// Define the GausND structure\n",
    "struct GausND {\n",
    "\n",
    "   TVectorD X;\n",
    "   TVectorD Mu;\n",
    "   TMatrixDSym CovMat;\n",
    "\n",
    "   GausND( int dim ) :\n",
    "      X(TVectorD(dim)),\n",
    "      Mu(TVectorD(dim)),\n",
    "      CovMat(TMatrixDSym(dim) )\n",
    "   {}\n",
    "\n",
    "   double operator() (double *x, double *p) {\n",
    "      // 4 parameters\n",
    "      int dim = X.GetNrows();\n",
    "      int k = 0;\n",
    "      for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; }\n",
    "      for (int i = 0; i<dim; ++i) {\n",
    "         CovMat(i,i) = p[k]*p[k];\n",
    "         k++;\n",
    "      }\n",
    "      for (int i = 0; i<dim; ++i) {\n",
    "         for (int j = i+1; j<dim; ++j) {\n",
    "            // p now are the correlations N(N-1)/2\n",
    "               CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j));\n",
    "               CovMat(j,i) = CovMat(i,j);\n",
    "               k++;\n",
    "         }\n",
    "      }\n",
    "      if (debug) {\n",
    "         X.Print();\n",
    "         CovMat.Print();\n",
    "      }\n",
    "\n",
    "      double det = CovMat.Determinant();\n",
    "      if (det <= 0) {\n",
    "         Fatal(\"GausND\",\"Determinant is <= 0 det = %f\",det);\n",
    "         CovMat.Print();\n",
    "         return 0;\n",
    "      }\n",
    "      double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det);\n",
    "      // compute the gaussians\n",
    "      CovMat.Invert();\n",
    "      double fval  = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm;\n",
    "\n",
    "      if (debug) {\n",
    "         std::cout << \"det  \" << det << std::endl;\n",
    "         std::cout << \"norm \" << norm << std::endl;\n",
    "         std::cout << \"fval \" << fval << std::endl;\n",
    "      }\n",
    "\n",
    "      return fval;\n",
    "   }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fb88c3dd",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:11.699833Z",
     "iopub.status.busy": "2026-05-19T20:26:11.699716Z",
     "iopub.status.idle": "2026-05-19T20:26:11.904681Z",
     "shell.execute_reply": "2026-05-19T20:26:11.903929Z"
    }
   },
   "outputs": [],
   "source": [
    "const int N = 10000;\n",
    "/*const int NBin = 1000;*/\n",
    "const int DIM = 4;\n",
    "\n",
    "double xmin[] = {-10,-10,-10, -10};\n",
    "double xmax[] = { 10, 10, 10,  10};\n",
    "double par0[] = { 1., -1., 2, 0, // the gaussian mu\n",
    "                  1, 2, 1, 3, // the sigma\n",
    "                  0.5,0.,0.,0.,0.,0.8 };  // the correlation\n",
    "\n",
    "const int NPAR = DIM + DIM*(DIM+1)/2; // 14 in the 4d case"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "823c461a",
   "metadata": {},
   "source": [
    "generate the sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5828910d",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:11.914451Z",
     "iopub.status.busy": "2026-05-19T20:26:11.914317Z",
     "iopub.status.idle": "2026-05-19T20:26:12.121998Z",
     "shell.execute_reply": "2026-05-19T20:26:12.121587Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_55:2:2: error: unknown type name 'GausND'\n",
      " GausND gaus4d(4);\n",
      " ^\n"
     ]
    }
   ],
   "source": [
    "GausND gaus4d(4);\n",
    "TF1 * f = new TF1(\"functionND\",gaus4d,0,1,14);\n",
    "f->SetParameters(par0);\n",
    "\n",
    "double x0[] = {0,0,0,0};"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1136beb",
   "metadata": {},
   "source": [
    "for debugging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f3c77dac",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:12.123987Z",
     "iopub.status.busy": "2026-05-19T20:26:12.123835Z",
     "iopub.status.idle": "2026-05-19T20:26:12.327740Z",
     "shell.execute_reply": "2026-05-19T20:26:12.327060Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_56:3:1: warning: declaration without the 'auto' keyword is deprecated: function '__cling_Un1Qu320' [-Wdeprecated-declarations]\n",
      "debug = false;\n",
      "^~~~~\n",
      "auto \n",
      "input_line_56:8:64: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const int *'\n",
      "   else if (i < 2*DIM) f->SetParName(i, name.Format(\"sig_%d\",i-DIM+1) );\n",
      "                                                               ^~~\n",
      "input_line_56:9:66: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const int *'\n",
      "   else if (i < 2*DIM) f->SetParName(i, name.Format(\"sig_%d\",i-2*DIM+1) );\n",
      "                                                                 ^~~\n",
      "input_line_56:25:25: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const int *'\n",
      "sampler->SetFunction(*f,DIM);\n",
      "                        ^~~\n"
     ]
    }
   ],
   "source": [
    "if (debug) f->EvalPar(x0,nullptr);\n",
    "debug = false;\n",
    "\n",
    "TString name;\n",
    "for (int i = 0; i < NPAR; ++i )  {\n",
    "   if (i < DIM) f->SetParName(i, name.Format(\"mu_%d\",i+1) );\n",
    "   else if (i < 2*DIM) f->SetParName(i, name.Format(\"sig_%d\",i-DIM+1) );\n",
    "   else if (i < 2*DIM) f->SetParName(i, name.Format(\"sig_%d\",i-2*DIM+1) );\n",
    "}\n",
    "\n",
    "/*ROOT::Math::DistSamplerOptions::SetDefaultSampler(\"Foam\");*/\n",
    "DistSampler * sampler = Factory::CreateDistSampler();\n",
    "if (sampler == nullptr) {\n",
    "   Info(\"multidimSampling\",\"Default sampler %s is not available try with Foam \",\n",
    "        ROOT::Math::DistSamplerOptions::DefaultSampler().c_str() );\n",
    "   ROOT::Math::DistSamplerOptions::SetDefaultSampler(\"Foam\");\n",
    "}\n",
    "sampler = Factory::CreateDistSampler();\n",
    "if (sampler == nullptr) {\n",
    "   Error(\"multidimSampling\",\"Foam sampler is not available - exit \");\n",
    "   return;\n",
    "}\n",
    "\n",
    "sampler->SetFunction(*f,DIM);\n",
    "sampler->SetRange(xmin,xmax);\n",
    "bool ret = sampler->Init();\n",
    "\n",
    "std::vector<double> data1(DIM*N);\n",
    "double v[DIM];\n",
    "TStopwatch w;\n",
    "\n",
    "if (!ret) {\n",
    "   Error(\"Sampler::Init\",\"Error initializing unuran sampler\");\n",
    "   return;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e59dd333",
   "metadata": {},
   "source": [
    "generate the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bf69fef6",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:12.329432Z",
     "iopub.status.busy": "2026-05-19T20:26:12.329313Z",
     "iopub.status.idle": "2026-05-19T20:26:12.536534Z",
     "shell.execute_reply": "2026-05-19T20:26:12.535831Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_57:6:13: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const int *'\n",
      "      data1[N*j + i]     = v[j];\n",
      "            ^\n"
     ]
    }
   ],
   "source": [
    "w.Start();\n",
    "for (int i = 0; i < N; ++i) {\n",
    "   sampler->Sample(v);\n",
    "   for (int j = 0; j < DIM; ++j)\n",
    "      data1[N*j + i]     = v[j];\n",
    "}\n",
    "w.Stop();\n",
    "w.Print();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08f5130d",
   "metadata": {},
   "source": [
    "fill tree with data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f51c748b",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:12.538301Z",
     "iopub.status.busy": "2026-05-19T20:26:12.538184Z",
     "iopub.status.idle": "2026-05-19T20:26:12.740773Z",
     "shell.execute_reply": "2026-05-19T20:26:12.740059Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_58:8:22: error: cannot initialize an array element of type 'void *' with an rvalue of type 'const int *'\n",
      "      x[j] = data1[i+N*j];\n",
      "                     ^\n"
     ]
    }
   ],
   "source": [
    "TFile * file = new TFile(\"multiDimSampling.root\",\"RECREATE\");\n",
    "double x[DIM];\n",
    "TTree * t1 = new TTree(\"t1\",\"Tree from Unuran\");\n",
    "t1->Branch(\"x\",x,\"x[4]/D\");\n",
    "for (int i = 0; i < N; ++i) {\n",
    "   for (int j = 0; j < DIM; ++j) {\n",
    "      x[j] = data1[i+N*j];\n",
    "   }\n",
    "   t1->Fill();\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6e5c499",
   "metadata": {},
   "source": [
    "plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "38101d1b",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:12.742552Z",
     "iopub.status.busy": "2026-05-19T20:26:12.742433Z",
     "iopub.status.idle": "2026-05-19T20:26:12.949780Z",
     "shell.execute_reply": "2026-05-19T20:26:12.949025Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_60:2:3: error: use of undeclared identifier 't1'\n",
      " (t1->Draw(\"x[0]:x[1]:x[2]:x[3]\", \"\", \"candle\"))\n",
      "  ^\n",
      "Error in <HandleInterpreterException>: Error evaluating expression (t1->Draw(\"x[0]:x[1]:x[2]:x[3]\", \"\", \"candle\"))\n",
      "Execution of your code was aborted.\n"
     ]
    }
   ],
   "source": [
    "t1->Draw(\"x[0]:x[1]:x[2]:x[3]\",\"\",\"candle\");\n",
    "TCanvas * c2 = new TCanvas();\n",
    "c2->Divide(3,2);\n",
    "int ic=1;\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[0]:x[1]\");\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[0]:x[2]\");\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[0]:x[3]\");\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[1]:x[2]\");\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[1]:x[3]\");\n",
    "c2->cd(ic++);\n",
    "t1->Draw(\"x[2]:x[3]\");\n",
    "\n",
    "t1->Write();\n",
    "file->Close();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d763140",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3e635f48",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2026-05-19T20:26:12.951503Z",
     "iopub.status.busy": "2026-05-19T20:26:12.951387Z",
     "iopub.status.idle": "2026-05-19T20:26:13.158502Z",
     "shell.execute_reply": "2026-05-19T20:26:13.157657Z"
    }
   },
   "outputs": [],
   "source": [
    "gROOT->GetListOfCanvases()->Draw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ROOT C++",
   "language": "c++",
   "name": "root"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".C",
   "mimetype": " text/x-c++src",
   "name": "c++"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}