FittingTest_TMinuit_Din_vz0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "08a39b26",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ROOT\n",
    "from ROOT import *\n",
    "import sys\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "42b9068d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_49:31:44: error: use of undeclared identifier 'my2Dfunc'\n",
      "         tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,(double*) p))/h1->GetBinError(ix,iy);\n",
      "                                           ^\n",
      "input_line_49:42:44: error: use of undeclared identifier 'my2Dfunc'\n",
      "         tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,(double *) p))/h2->GetBinError(ix,iy);\n",
      "                                           ^\n"
     ]
    }
   ],
   "source": [
    "%%cpp -d\n",
    "class MyFcn {\n",
    " public:\n",
    "   TH2D *h1 = nullptr;\n",
    "   TH2D *h2 = nullptr;\n",
    "   int npfits = 0;\n",
    "\n",
    "   MyFcn(TH2D * _h1, TH2D * _h2) :\n",
    "    h1(_h1), h2(_h2) {}\n",
    "\n",
    "   double operator()(const double *p) {\n",
    "\n",
    "   TAxis *xaxis1  = h1->GetXaxis();\n",
    "   TAxis *yaxis1  = h1->GetYaxis();\n",
    "   TAxis *xaxis2  = h2->GetXaxis();\n",
    "   TAxis *yaxis2  = h2->GetYaxis();\n",
    "\n",
    "   int nbinX1 = h1->GetNbinsX();\n",
    "   int nbinY1 = h1->GetNbinsY();\n",
    "   int nbinX2 = h2->GetNbinsX();\n",
    "   int nbinY2 = h2->GetNbinsY();\n",
    "\n",
    "   double chi2 = 0;\n",
    "   double x[2];\n",
    "   double tmp;\n",
    "   npfits = 0;\n",
    "   for (int ix = 1; ix <= nbinX1; ++ix) {\n",
    "      x[0] = xaxis1->GetBinCenter(ix);\n",
    "      for (int iy = 1; iy <= nbinY1; ++iy) {\n",
    "         if ( h1->GetBinError(ix,iy) > 0 ) {\n",
    "         x[1] = yaxis1->GetBinCenter(iy);\n",
    "         tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,(double*) p))/h1->GetBinError(ix,iy);\n",
    "         chi2 += tmp*tmp;\n",
    "         npfits++;\n",
    "         }\n",
    "      }\n",
    "   }\n",
    "   for (int ix = 1; ix <= nbinX2; ++ix) {\n",
    "      x[0] = xaxis2->GetBinCenter(ix);\n",
    "      for (int iy = 1; iy <= nbinY2; ++iy) {\n",
    "         if ( h2->GetBinError(ix,iy) > 0 ) {\n",
    "         x[1] = yaxis2->GetBinCenter(iy);\n",
    "         tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,(double *) p))/h2->GetBinError(ix,iy);\n",
    "         chi2 += tmp*tmp;\n",
    "         npfits++;\n",
    "         }\n",
    "      }\n",
    "   }\n",
    "   return chi2;\n",
    "}\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ac1a82a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "#include \"TH2D.h\"\n",
    "#include \"TF2.h\"\n",
    "#include \"TCanvas.h\"\n",
    "#include \"TStyle.h\"\n",
    "#include \"TRandom3.h\"\n",
    "#include \"Fit/Fitter.h\"\n",
    "#include \"TList.h\"\n",
    "\n",
    "#include <iostream>\n",
    "\n",
    "double gauss2D(double *x, double *par) {\n",
    "   double z1 = double((x[0]-par[1])/par[2]);\n",
    "   double z2 = double((x[1]-par[3])/par[4]);\n",
    "   return par[0]*exp(-0.5*(z1*z1+z2*z2));\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "94d9a9aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "double my2Dfunc(double *x, double *par) {\n",
    "   return gauss2D(x,&par[0]) + gauss2D(x,&par[5]);\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "498d2433",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "void FillHisto(TH2D * h, int n, double * p) {\n",
    "\n",
    "\n",
    "   const double mx1 = p[1];\n",
    "   const double my1 = p[3];\n",
    "   const double sx1 = p[2];\n",
    "   const double sy1 = p[4];\n",
    "   const double mx2 = p[6];\n",
    "   const double my2 = p[8];\n",
    "   const double sx2 = p[7];\n",
    "   const double sy2 = p[9];\n",
    "   //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2);\n",
    "   const double w1 = 0.5;\n",
    "\n",
    "   double x, y;\n",
    "   for (int i = 0; i < n; ++i) {\n",
    "      // generate randoms with larger Gaussians\n",
    "      gRandom->Rannor(x,y);\n",
    "\n",
    "      double r = gRandom->Rndm(1);\n",
    "      if (r < w1) {\n",
    "         x = x*sx1 + mx1;\n",
    "         y = y*sy1 + my1;\n",
    "      }\n",
    "      else {\n",
    "         x = x*sx2 + mx2;\n",
    "         y = y*sy2 + my2;\n",
    "      }\n",
    "      h->Fill(x,y);\n",
    "\n",
    "   }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "813e985e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "int option=1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "14498be1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "int nbx1 = 50;\n",
    "int nby1 = 50;\n",
    "int nbx2 = 50;\n",
    "int nby2 = 50;\n",
    "double xlow1 = 0.;\n",
    "double ylow1 = 0.;\n",
    "double xup1 = 10.;\n",
    "double yup1 = 10.;\n",
    "double xlow2 = 5.;\n",
    "double ylow2 = 5.;\n",
    "double xup2 = 20.;\n",
    "double yup2 = 20.;\n",
    "\n",
    "auto h1 = new TH2D(\"h1\",\"core\",nbx1,xlow1,xup1,nby1,ylow1,yup1);\n",
    "auto h2 = new TH2D(\"h2\",\"tails\",nbx2,xlow2,xup2,nby2,ylow2,yup2);\n",
    "\n",
    "double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8114216e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_55:3:1: error: unknown type name 'func'\n",
      "func->SetParameters(iniParams);\n",
      "^\n",
      "input_line_55:3:5: error: cannot use arrow operator on a type\n",
      "func->SetParameters(iniParams);\n",
      "    ^\n"
     ]
    }
   ],
   "source": [
    "%%cpp -d\n",
    "\n",
    "TF2 * func = new TF2(\"func\",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);\n",
    "func->SetParameters(iniParams);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "135ef9b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "input_line_56:4:1: error: C++ requires a type specifier for all declarations\n",
      "FillHisto(h1,n1,iniParams);\n",
      "^\n",
      "input_line_56:5:1: error: C++ requires a type specifier for all declarations\n",
      "FillHisto(h2,n2,iniParams);\n",
      "^\n"
     ]
    }
   ],
   "source": [
    "%%cpp -d\n",
    "\n",
    "int n1 = 1000000;\n",
    "int n2 = 1000000;\n",
    "FillHisto(h1,n1,iniParams);\n",
    "FillHisto(h2,n2,iniParams);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92a226b5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}