Newer
Older
R_phipi / minuit_fitter.ipynb
@Davide Lancierini Davide Lancierini on 3 Oct 2018 4 KB added minuit fitter script
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ROOT as r\n",
    "import ctypes\n",
    "import numpy as np\n",
    "from array import array\n",
    "import root_numpy as rn\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fcn( npar, gin, f, par, iflag):\n",
    "    \n",
    "    global ncount\n",
    "    global chisq\n",
    "    global pulls\n",
    "    global fitfunavg\n",
    "    global delta\n",
    "    \n",
    "    pull = np.array([0. for i in range(nbins)])\n",
    "    residuals = np.array([0. for i in range(nbins)])\n",
    "    \n",
    "    fitfunavg=np.array([func((xup[i]+xdown[i])/2) for i in range(nbins)])\n",
    "    \n",
    "    chisq, delta = 0., 0.\n",
    "    \n",
    "    delta = np.array([data[i]-fitfun[i] for i in range(0,nbins)])\n",
    "    \n",
    "    chisq = np.dot(np.transpose(delta),np.dot(Vxinv,delta))\n",
    "    \n",
    "    for i in range(nbins):\n",
    "        \n",
    "        pulls[i] = (y[i]-fitfun[i])/(np.sqrt(Vx[i,i]))\n",
    "        \n",
    "    print('X^2', chisq)\n",
    "    f[0] = chisq\n",
    "    ncount+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###DEFINE FIT FUNCTIONS HERE\n",
    "###AS FUNC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##___MINUIT______\n",
    "\n",
    "def fit():\n",
    "    \n",
    "    gMinuit = r.TMinuit(5)\n",
    "    gMinuit.SetFCN(fcn)\n",
    "    \n",
    "    arglist = np.array([0. for i in range(10)], dtype='float32')\n",
    "    ierflg = ctypes.c_int(1982)\n",
    "    \n",
    "    arglist[0]=1\n",
    "    gMinuit.mnexcm( \"SET ERR\", arglist, 1, ierflg )\n",
    "    \n",
    "    # Set starting values and step sizes for parameters\n",
    "    vstart = np.array([], dtype='float32')\n",
    "    step =  np.array([], dtype='float32')\n",
    "    \n",
    "    gMinuit.mnparm( 0, \"\", vstart[0], step[0], lower_bound, upper_bound, ierflg)#par0\n",
    "    \n",
    "    #if needed\n",
    "    #gMinuit.FixParameter(0)\n",
    "    \n",
    "    #Minimization Step\n",
    "    \n",
    "    arglist[0] = 500\n",
    "    arglist[1] = 1e-20\n",
    "    gMinuit.mnexcm( \"MIGRAD\", arglist, 3, ierflg)\n",
    "    gMinuit.mnexcm(\"HESSE\", arglist, 4, ierflg)\n",
    "    \n",
    "    amin, edm, errdef, par0, err_par0 = map(ctypes.c_double, (0.18, 0.19, 0.20, 0.21, 0.22, 0.23))\n",
    "    nvpar, nparx, icstat = map(ctypes.c_int, (1983, 1984, 1985))\n",
    "    gMinuit.mnstat( amin, edm, errdef, nvpar, nparx, icstat )\n",
    "    gMinuit.mnprin( 4, amin.value )\n",
    "    gMinuit.GetParameter(0,par0,err_par0)\n",
    "    \n",
    "    fitpars = np.array([par0.value])\n",
    "    err_fitpars = np.array([err_par0.value])\n",
    "    ratio = chisq/(fit_bins-2)\n",
    "    \n",
    "    return firpars, err_fitpars, ratio, pulls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__maun__':\n",
    "    \n",
    "    binwidth = (end-start)/nbins\n",
    "    \n",
    "    xup = np.array([start +1.0*(i+1.)*binwidth] for i in range(0,nbins))\n",
    "    xdown = np.array([start +1.0*(i+0.)*binwidth] for i in range(0,nbins))\n",
    "    \n",
    "    fitpars, err_fitpars, ratio, pulls = fit()\n",
    "    \n",
    "    #save data and plot here\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}