{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Libraries needed by the following routines.\n", "import numpy as np\n", "from iminuit import Minuit\n", "# Functions we'll need that must be imported from external libraries.\n", "from math import sqrt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Global variables that are used for inter-routine communication. Note that \"global\" just makes the variables accessible elsewhere in the Python program. It does not define values for them.\n", "\n", "By the way, globals are considered poor practice in Python programming. However, in translating the FORTRAN code and COMMON blocks, it's the simplest solution for now." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Define global variables\n", "global dat_array\n", "global sinsqthw0\n", "global nev\n", "global evt_evis \n", "global evt_evis0\n", "global bkg_evis0\n", "global bkg_IBD_evis0\n", "global bkg_total\n", "global bkg_IBD_total\n", "global evis_min\n", "global evis_max\n", "global d_evis\n", "global n_evis\n", "global i_evis\n", "global istat_only\n", "\n", "global f_res\n", "global gsig\n", "global ngsig\n", "\n", "global eps_eL_ee0\n", "global eps_eR_ee0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A reminder of the COMMON blocks in the original FORTRAN program:\n", "\n", "```\n", "# data array: enu,levt,e_evis,wtnue,wtnuebar\n", "real*8 sinsqthw0\n", "integer nev,nev_max/2000000/\n", "real*8 dat_array(5,2000000)\n", "common /data/dat_array,sinsqthw0,nev\n", "\n", "# Evis bins and limits\n", "integer n_evis,i_evis,istat_only\n", "real*8 evis_min,evis_max,d_evis\n", "real*8 evt_evis(100),evt_evis0(100)\n", "real*8 bkg_evis0(100),bkg_total\n", "real*8 bkg_IBD_evis0(100),bkg_IBD_total\n", "common /evisbins/evt_evis,evt_evis0,bkg_evis0,\n", "> bkg_IBD_evis0,bkg_total,bkg_IBD_total,\n", "> evis_min,evis_max,d_evis,n_evis,i_evis,istat_only\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function that we'll fit using Minuit." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def fcn(sinsqthw,norm,normIBD,normBkgnd,eps_eR_ee,eps_eL_ee,eps_eR_te,eps_eL_te,evis_res,gL2,gR2,gL,gR):\n", " \n", " global sinsqthw0\n", " global eps_eL_ee0\n", " global eps_eR_ee0\n", " global dat_array\n", " global istat_only\n", " global evis_min\n", " global evis_max\n", " global d_evis\n", " global n_evis\n", " global i_evis\n", " global gsig\n", " global ngsig\n", " global f_res\n", " global bkg_evis0\n", " global bkg_IBD_evis0\n", " global evt_evis0\n", " \n", " d_norm = 0.007 # <====== Nominal = 0.1% or 0.001 <== IBD eff uncertainty = 0.7%\n", " d_normIBD = 0.08 # <===== New eff = 0.9975+/-0.0002 => dnorm = 0.0002/0.0025 = 0.08\n", " d_normBkgnd = 1.0/sqrt(bkg_total) # <== Use beam off measurements but use bin by bin errors##\n", " d_evis_res = 0.001 # evis energy resolution uncertainty 6.5% +/- 0.1%\n", " d_eps_eR_ee = 3.0\n", " d_eps_eL_ee = 3.0\n", " d_eps_eR_te = 3.0\n", " d_eps_eL_te = 3.0\n", "\n", " # shape only\n", " # d_norm = 10.\n", " \n", " global first_fcn\n", " try:\n", " first_fcn\n", " except NameError:\n", " print (\"d_norm,d_normIBD: {:7.4f}{:7.4f}\".format(d_norm,d_normIBD))\n", " first_fcn = 1\n", "\n", " # add pull terms if required\n", " chi2 = (((1.-norm)/d_norm)**2 \n", " + ((1.-normBkgnd)/d_normBkgnd)**2\n", " + ((1.-normIBD)/d_normIBD)**2\n", " # + (eps_eR_ee/d_eps_eR_ee)**2\n", " # + (eps_eL_ee/d_eps_eL_ee)**2\n", " # + (eps_eR_te/d_eps_eR_te)**2\n", " # + (eps_eL_te/d_eps_eL_te)**2\n", " + (evis_res/d_evis_res)**2)\n", "\n", "\n", " # set up helicity factors including possible NSI factors\n", " eminus0 = -0.5-sinsqthw0-eps_eL_ee0\n", " eplus0 = -sinsqthw0-eps_eR_ee0\n", "\n", " eminus = -0.5-sinsqthw-eps_eL_ee\n", " eplus = -sinsqthw-eps_eR_ee\n", "\n", " # option to fit for gL2 and gR2\n", " # eminus = -sqrt(gL2)/2.\n", " # eplus = -sqrt(gR2)/2.\n", "\n", " # option to fit for gL and gR\n", " eminus = -gL\n", " eplus = -gR\n", "\n", " # bin data with new sinsqthw\n", " # data_array: enu,levt,e_evis,wtnue,wtnuebar\n", " # zero out evt_evis array\n", " evt_evis = np.zeros((n_evis))\n", " # reweight events for new sinsqthw and bin events in evis\n", " for iev in range(0,nev):\n", " evis = dat_array[3,iev]\n", " enu = dat_array[1,iev]\n", " y = evis/enu\n", " factor0 = (eplus0**2+eminus0**2*(1-y)**2\n", " - eminus0*eplus0*y*0.511/enu) # This term changes 0.008425 to 0.07960 <== Check sign\n", " # / (eplus0**2+eminus0**2/3.0) <== Mistake in reweighting\n", " factor =( (eplus**2+eps_eR_te**2)\n", " + (eminus**2+eps_eL_te**2)*(1-y)**2 \n", " - eminus*eplus*y*0.511/enu) # This term changes 0.008425 to 0.07960 <== Check sign\n", " # / ((eplus**2+eps_eR_te**2) <== Mistake in reweighting\n", " # +(eminus**2+eps_eL_te**2)/3.0)\n", " for ig in range(0,ngsig):\n", " wtnuebar = dat_array[4,iev]/float(ngsig)\n", " esig = (f_res+evis_res)*sqrt(dat_array[3,iev])\n", " evis = dat_array[2,iev] + gsig[ig]*esig\n", " i_evis = int((evis-evis_min)/d_evis)\n", " i_evis = min(n_evis-1,max(0,i_evis))\n", " if (evis >= evis_min and evis <= evis_max):\n", " evt_evis[i_evis] += wtnuebar*(factor/factor0)\n", " \n", " # calculate chi2 over visible energy bins\n", " for i_evis in range(0,n_evis):\n", " pred = ( norm*evt_evis[i_evis]\n", " + normBkgnd*bkg_evis0[i_evis]\n", " + normIBD*bkg_IBD_evis0[i_evis] )\n", " data = ( evt_evis0[i_evis]\n", " + bkg_evis0[i_evis] + bkg_IBD_evis0[i_evis] )\n", " if (data > 10.):\n", " if istat_only:\n", " chi2 += (data-pred)**2/max(1.,data)\n", " else:\n", " chi2 += ( (data-pred)**2 \n", " /max(1.,data+bkg_evis0[i_evis]) ) # <== syst bkgnd error = stat error bkgnd\n", "\n", " # use binned max likelihood statistic\n", " # 2*(pred-data-data*log(pred/data))\n", " \n", " return chi2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These variables were defined identically in the main program and in fcn. I've taken the liberty of defining them only once, and making them global. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# energy smearing\n", "f_res = 0.065\n", "gsig = [-1.4,-0.53,0.0,0.53,1.4]\n", "# gsig = [0.0] # no smearing\n", "ngsig = len(gsig)\n", "\n", "eps_eL_ee0 = 0.0\n", "eps_eR_ee0 = 0.0" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "evis_min,evis_max,d_evis,n_evis: 3.0 12.0 0.2 45\n", "10000 1000000 0.238\n", "\n", "total flux: 1.147e+23 nyrs(with DF): 4.00 protons/yr: 1.97e+24\n", "\n", "isotopes/proton: 0.01456 isotope_br: 1.00\n", "\n", "beta-decay endpoint = 15.000 MeV\n", "\n", "mass: 897.0 radius: 6.50 shield: 3.2 distance: 12.50\n", "\n" ] } ], "source": [ "rad_run = [6.0,5.5,5.0]\n", "nrun = len(rad_run)\n", "\n", "# setup\n", "istat_only = False # =True for stat_only , =False for stat + syst\n", "twopar = True # Flag for two parameter fit (change dchi2 = 1.0 => 2.296)\n", "bkg_reduct = 1.0 # <===== Option bkg reduction by directionality\n", "radius0 = 6.5\n", "df_factor = 0.9/0.8 # 80%->90% DF for ES only (bkgnds already have this)\n", "veto_livetime = 0.624 # veto livetime\n", "irun = 2 # <== Standardize on radius<5m only (irun=2)\n", "radius = rad_run[irun]\n", "evis_min = 3.0 # MeV lower evis cut ## Also see below\n", "evis_max = 12.0 # MeV high evis cut (nominal 11 MeV) <== Matt wants 7 MeV for bkgnd\n", "d_evis = 0.2 # MeV evis bin size\n", "n_evis = int((evis_max-evis_min)/d_evis)\n", "\n", "print 'evis_min,evis_max,d_evis,n_evis: ', \\\n", " evis_min,evis_max,d_evis,n_evis\n", " \n", "# print out run parameters\n", "with open(\"current_istope_event_calc_run.txt\",\"r\") as f:\n", " run_parameters = f.readlines() # Read the file into a list of lines.\n", "for i in range(0,5):\n", " print run_parameters[i]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Python syntax note: For these multi-line equations, Python will not interpret them correctly unless there's a parenthesis in the first line to indicate there must be a continuation. For example, this will not work:\n", "\n", "```\n", "a = b\n", " + c\n", " + d\n", " + e\n", "```\n", "\n", "but this will:\n", "\n", "```\n", "a = ( b\n", " + c\n", " + d\n", " + e )\n", "```\n", "\n", "You can also use a continuation marker (\\\\) at the end of a line, but you have to make there is *no* space following the \\:\n", "\n", "```\n", "a = b \\\n", " + c \\\n", " + d \\\n", " + e\n", "```" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nev, sinsqthw0 = 1000000 0.238\n" ] } ], "source": [ "# read in elastic scattering data array\n", "with open(\"nuebar_evis_array.dat\",\"r\") as f:\n", " contents = f.readlines() # Read the file into a list of lines.\n", "# The first line contains the number of lines in the rest of the file.\n", "line = contents[0].strip()\n", "fields = line.split()\n", "# There's a mis-match between the actual number of lines in this file\n", "# and the number in the file, so get 'nev' from the actual number of\n", "# entires in the file (omitting the first line).\n", "nev = len(contents)-1 \n", "# sinsqthw0 is the second non-blank field. \n", "sinsqthw0 = float(fields[1])\n", "print 'nev, sinsqthw0 = ',nev,\" \",sinsqthw0\n", "\n", "# option to change the input distribution\n", "eminus0 = -0.5-sinsqthw0\n", "eplus0 = -sinsqthw0\n", "eminus = -0.5-sinsqthw0-eps_eL_ee0\n", "eplus = -sinsqthw0-eps_eR_ee0\n", "dat_array = np.zeros((5,nev))\n", "for i in range(1,nev+1):\n", " iev = i-1 # Adjust for the first line of the file\n", " # dat_array[k,iev] where k = enu,levt,e_evis,wtnue,wtnuebar\n", " fields = contents[i].split()\n", " for j in range(0,5):\n", " dat_array[j,iev] = float(fields[j])\n", "\n", " evis = dat_array[2,iev]\n", " enu = dat_array[0,iev]\n", " y = evis/enu\n", " factor0 = eplus0**2 + eminus0**2*(1-y)**2\n", " factor = eplus**2 + eminus**2*(1-y)**2 \\\n", " + eminus*eplus*y*0.511/enu # This term changes 0.008425 to 0.07960\n", " # rescale for radius\n", " dat_array[3,iev] *= veto_livetime*df_factor*(radius/radius0)**3 \\\n", " * factor/factor0\n", " dat_array[4,iev] *= veto_livetime*df_factor*(radius/radius0)**3 \\\n", " * factor/factor0" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "IBD_2.txt\n" ] } ], "source": [ "# IBD background is separate since tied to cyclotron\n", "IBD_name = 'IBD_' # IBD background is special since tied to cyclotron\n", "nebk = 160\n", "bkg_IBD = np.zeros((nebk))\n", "evis_IBD = np.zeros((nebk))\n", "\n", "suffix = str(irun)\n", "filename = IBD_name + suffix + '.txt'\n", "print filename\n", "with open(filename,\"r\") as f:\n", " contents = f.readlines() # Read the file into a list of lines.\n", "for iebk in range(0, nebk):\n", " fields = contents[iebk].split()\n", " e = float(fields[0])\n", " evt = float(fields[1])\n", " evt *= bkg_reduct # <===== Option bkg reduction by directionality\n", " # evt *= 0.25 # <= New IBD eff = 0.9975+/-0.0002 => Old was 0.99 <= New IBD file Toups not needed 7/12/13\n", " # evt *= df_factor # <===== Correct for new duty factor <= Files already have 90% df 7/12/13\n", " bkg_IBD[iebk] += evt\n", " evis_IBD[iebk] = e" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "B8Solar_2.txt\n", "8Li_2.txt\n", "8B_2.txt\n", "11Be_2.txt\n", "gamma_stain_2.txt\n", "gamma_rock_2.txt\n", "208Tl_2.txt\n", " evis bk_total B8Solar_ 8Li_ 8B_ 11Be_ gamma_stain_ gamma_rock_ 208Tl_ IBD_\n", " 0.05 31.56 30.48 0 0 1.08 0 0 0 0\n", " 0.15 32.25 30.26 0 0 1.99 0 0 0 0\n", " 0.25 32.78 30.05 0 0 2.735 0 0 0 0\n", " 0.35 33.27 29.85 0.01783 0 3.403 0 0 0 0\n", " 0.45 33.7 29.65 0.05311 0 3.996 0 0 0 0\n", " 0.55 34.04 29.45 0.0801 0 4.511 0 0 0 0\n", " 0.65 34.31 29.26 0.1083 0 4.945 0 0 0 0\n" ] } ], "source": [ "# read in other backgrounds\n", "# nominal is 7 bkgnd types\n", "bkg_name =['B8Solar_','8Li_','8B_','11Be_','gamma_stain_','gamma_rock_','208Tl_']\n", "nbkg = len(bkg_name)\n", "evis_bkg = np.zeros((nebk))\n", "bkg_tot = np.zeros((nebk))\n", "bkg_types = np.zeros((nebk,nbkg))\n", "\n", "for ibkg in range(0,nbkg):\n", " filename = bkg_name[ibkg] + suffix + \".txt\"\n", " print filename\n", " with open(filename,\"r\") as f:\n", " contents = f.readlines() # Read the file into a list of lines.\n", " for iebk in range(0,nebk):\n", " fields = contents[iebk].split()\n", " e = float(fields[0])\n", " evt = float(fields[1])\n", " evt *= bkg_reduct # <===== Option bkg reduction by directionality\n", " # evt *= df_factor # <===== Correct for new duty factor <= Files already have 90% df 7/12/13\n", " bkg_types[iebk,ibkg] += evt\n", " bkg_tot[iebk] += evt\n", " evis_bkg[iebk] = e # <== this assumes that all backgnds have same evis binning\n", "\n", "# In Python, if your print statement ends in a comma, it won't\n", "# append a end-of-line on output.\n", "# {:>12s} means to print the string right-aligned within the 12-character field.\n", "print \"{:>12s}{:>12s}\".format(\"evis\",\"bk_total\"),\n", "for i in range(0,nbkg):\n", " print \"{:>12s}\".format(bkg_name[i]),\n", "print \"IBD_\"\n", "\n", "for iebk in range(0,nbkg):\n", " print \"{:12.4g}{:12.4g}\".format(evis_bkg[iebk],bkg_tot[iebk]),\n", " for ibkg in range (0,nbkg):\n", " print \"{:12.4g}\".format(bkg_types[iebk,ibkg]),\n", " print \"{:12.4g}\".format(bkg_IBD[iebk])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the original FORTRAN, this was given as:\n", "\n", "```\n", "do icut=3,3\n", "```\n", "\n", "For now, since I'm developing the Python code section-by-section, I'm just setting icode to 3. When moving to cut studies, the following sections will have to put into a Python 'for' loop. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "icut = 3" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "evis_min,evis_max,evt_total: 3.0 12.0 4222.18363683\n", "evt_total,bkg_total,bk_IBD_total: 4222.18363683 2870.0131834 705.261233\n", "evis,evt_evis0,bkg_evis0,bkg_IBD_evis0\n", "3.1 352.600859296 105.5690582 5.7753656\n", "3.3 327.170537321 155.3371452 8.8630442\n", "3.5 300.88692453 193.0240474 9.8454106\n", "3.7 280.548308795 191.183455 10.7459636\n", "3.9 258.847331367 175.0627758 11.8161394\n", "4.1 238.766072412 151.0194352 12.881082\n", "4.3 218.770612887 124.8238746 13.8414446\n", "4.5 204.549367005 103.3272248 14.7774606\n", "4.7 187.472602224 91.8323788 15.750869\n", "4.9 172.203007557 87.1322806 16.7545912\n", "5.1 159.307558471 83.916173 17.7079412\n", "5.3 148.607800488 77.054905 18.5467566\n", "5.5 135.10437807 77.2525484 19.2209426\n", "5.7 122.793922875 76.5204998 19.9506278\n", "5.9 112.949476001 77.0357898 21.046389\n", "6.1 104.287030511 78.8436218 21.7312504\n", "6.3 95.5424534113 84.2867754 22.1445562\n", "6.5 87.8203520606 82.1292952 22.2499084\n", "6.7 77.7732122149 69.3695554 22.603908\n", "6.9 73.2538316653 67.6956312 23.0646532\n", "7.1 64.7789755504 68.2705032 23.3513514\n", "7.3 60.2845732547 86.7551774 25.4591692\n", "7.5 54.7233540091 75.9251138 23.156071\n", "7.7 48.6664341943 89.239677 20.6261012\n", "7.9 43.6831854347 64.8001656 22.7188408\n", "8.1 38.1899067707 42.7635478 22.2986882\n", "8.3 35.8844780942 34.3146172 21.6527688\n", "8.5 31.3291656298 37.6880908 21.211009\n", "8.7 27.9615501467 33.4713254 20.4867074\n", "8.9 24.6794952764 37.1682804 19.6691426\n", "9.1 21.6652772574 26.8593116 18.949852\n", "9.3 18.6094436308 23.659712 17.921198\n", "9.5 16.2825907851 18.1620792 16.907915\n", "9.7 13.8636370735 15.5038058 15.5409504\n", "9.9 12.5544579955 12.2019116 14.2721388\n", "10.1 10.3601029753 10.6917036 13.114181\n", "10.3 8.46698561059 8.2511002 11.9209164\n", "10.5 7.49408620029 7.0587494 10.5929004\n", "10.7 5.83398462404 6.4529204 8.976007\n", "10.9 4.96419102117 5.0537988 7.7055794\n", "11.1 3.90700761403 3.7703168 6.2403234\n", "11.3 3.09267938399 3.1059712 5.1449746\n", "11.5 2.4997219151 2.563913 3.813021\n", "11.7 1.80375498579 2.119919 2.612689\n", "11.9 1.34895822877 1.7450016 1.6004328\n" ] } ], "source": [ "# setup\n", "evis_min = float(icut) # MeV lower evis cut\n", "n_evis = int((evis_max-evis_min)/d_evis)\n", "\n", "# zero out evt_evis array\n", "evt_total = 0.\n", "bkg_total = 0.\n", "bkg_IBD_total = 0.\n", "evt_evis0 = np.zeros((n_evis))\n", "bkg_evis0 = np.zeros((n_evis))\n", "bkg_IBD_evis0 = np.zeros((n_evis))\n", "\n", "# bin ES events in evis\n", "for iev in range(0,nev):\n", " for ig in range(0,ngsig):\n", " enu = dat_array[0,iev]\n", " esig = f_res*sqrt(dat_array[2,iev])\n", " evis = dat_array[2,iev] + gsig[ig]*esig\n", " wtnuebar = dat_array[4,iev]/ngsig\n", " i_evis = int((evis-evis_min)/d_evis)\n", " i_evis = min(n_evis-1,max(0,i_evis))\n", " if (evis >= evis_min and evis <= evis_max):\n", " evt_evis0[i_evis] += wtnuebar\n", " evt_total += wtnuebar\n", "\n", "# bin non-beam bkg events in evis\n", "for iebk in range(0,nebk):\n", " for ig in range(0,ngsig):\n", " esig = f_res*sqrt(evis_bkg[iebk])\n", " evis = evis_bkg[iebk] + gsig[ig]*esig\n", " wtnuebar = bkg_tot[iebk]/ngsig\n", " i_evis = int((evis-evis_min)/d_evis)\n", " i_evis = min(n_evis-1,max(0,i_evis))\n", " if (evis >= evis_min and evis <= evis_max):\n", " bkg_evis0[i_evis] += wtnuebar\n", " bkg_total += wtnuebar\n", "\n", "# bin IBD (beam) bkg events in evis\n", "for iebk in range(0,nebk):\n", " for ig in range(0,ngsig):\n", " esig = f_res*sqrt(evis_IBD[iebk])\n", " evis = evis_IBD[iebk] + gsig[ig]*esig\n", " wtIBD = bkg_IBD[iebk]/ngsig\n", " i_evis = int((evis-evis_min)/d_evis)\n", " i_evis = min(n_evis-1,max(0,i_evis))\n", " if (evis >= evis_min and evis <= evis_max):\n", " bkg_IBD_evis0[i_evis] += wtIBD\n", " bkg_IBD_total += wtIBD\n", "\n", "if icut == 3:\n", " print 'evis_min,evis_max,evt_total: ',evis_min,evis_max,evt_total\n", " print 'evt_total,bkg_total,bk_IBD_total: ',evt_total,bkg_total,bkg_IBD_total\n", " print 'evis,evt_evis0,bkg_evis0,bkg_IBD_evis0'\n", " for i_evis in range(0,n_evis):\n", " evis = evis_min + (i_evis+0.5)*d_evis\n", " print evis, evt_evis0[i_evis],bkg_evis0[i_evis],bkg_IBD_evis0[i_evis]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The original FORTRAN program wrote to UNIT=12, and relied on FORTRAN's defaults to create a file name. Here I use `sinsqthw_nuebar_e_unit12.txt`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "unit12 = open(\"sinsqthw_nuebar_e_unit12.txt\",\"w\")\n", "\n", "# print out events vs sinsqthw\n", "unit12.write(\"sinsqthw,ES_evt_total\\n\")\n", "for isinsq in range(0,16):\n", " sinsqthw = sinsqthw0 + float(isinsq-7)*0.005\n", " eminus0 = -0.5-sinsqthw0\n", " eplus0 = -sinsqthw0\n", " eminus = -0.5-sinsqthw\n", " eplus = -sinsqthw\n", " evt_total_isinsq = 0.\n", " evt_evis = np.zeros((n_evis))\n", " # reweight events for new sinsqthw and bin events in evis\n", " for iev in range(0,nev):\n", " evis = dat_array[2,iev]\n", " enu = dat_array[0,iev]\n", " y = evis/enu\n", " factor0 = (eplus0**2+eminus0**2*(1-y)**2)\n", " # / (eplus0**2+eminus0**2/3.0)\n", " factor =( (eplus**2)\n", " + (eminus**2)*(1-y)**2)\n", " # / ((eplus**2)\n", " # +(eminus**2)/3.0)\n", " for ig in range(0,ngsig):\n", " wtnuebar = dat_array[4,iev]/ngsig\n", " esig = f_res*sqrt(dat_array[2,iev])\n", " evis = dat_array[2,iev] + gsig[ig]*esig\n", " i_evis = int((evis-evis_min)/d_evis)\n", " i_evis = min(n_evis-1,max(0,i_evis))\n", " if (evis >= evis_min and evis <= evis_max):\n", " evt_total_isinsq += wtnuebar*(factor/factor0)\n", " evt_evis[i_evis] += wtnuebar*(factor/factor0)\n", " unit12.write(str(sinsqthw) + \" \" + str(evt_total_isinsq) + \"\\n\")\n", " for i_evis in range(0,n_evis):\n", " evis = evis_min + (float(i_evis)+0.5)*d_evis\n", " # unit12.write(str(evis + \" \" + str(evt_evis[i_evis]) + \" \" + \\\n", " # + \" \" + str(bkg_evis0[i_evis]) + \" \" + str(bkg_IBD_evis0[i_evis]) + \"\\n\")\n", "\n", "unit12.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Python, we fiddle with the operation of Minuit by passing keyword arguments to the Minuit method. Because there are so many arguments for this code, we're using a Python dictionary to set the values. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FCN = 0.110482601489TOTAL NCALL = 54NCALLS = 54
EDM = 1.60226263358e-05GOAL EDM = 2.296e-05\n", " UP = 2.296
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ValidValid ParamAccurate CovarPosDefMade PosDef
TrueTrueFalseFalseTrue
Hesse FailHasCovAbove EDMReach calllim
FalseTrueFalseFalse
\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
+NameValueParab ErrorMinos Error-Minos Error+Limit-Limit+FIXED
1sinsqthw0.238100FIXED
2norm1.000020.010592800
3normIBD1.000510.10264900
4normBkgnd1100FIXED
5eps_eR_ee0100FIXED
6eps_eL_ee0100FIXED
7eps_eR_te0100FIXED
8eps_eL_te0100FIXED
9evis_res0100FIXED
10gL22.17857100FIXED
11gR20.226576100FIXED
12gL0.7381130.21499700
13gR0.2383460.65634400
\n", " \n", " 
\n",
       "            \n",
       "
\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# ============================================================================\n", "\n", "# initialize minuit\n", "# set the print level\n", "# print_level = 0 => Don't print anything\n", "# print_level = 1 => Print results at the end\n", "# print_level = 2 => Print fit status during the call\n", "arguments = dict(print_level=1)\n", "\n", "# set value for error calculation\n", "errdef = 1.0\n", "# Flag for two parameter fit (change dchi2 = 1.0 => 2.296)\n", "if twopar:\n", " errdef = 2.296\n", "arguments.update(dict(errordef=errdef))\n", "\n", "# set up parameters that vary\n", "arguments.update(dict(sinsqthw=sinsqthw0, fix_sinsqthw=True))\n", "if istat_only:\n", " arguments.update(dict(norm=1.0, fix_norm=True))\n", " arguments.update(dict(normIBD=1.0, fix_normIBD=True))\n", "else:\n", " arguments.update(dict(norm=1.0, error_norm=0.1))\n", " arguments.update(dict(normIBD=1.0, error_normIBD=0.1))\n", "\n", "# Put Bkgnd error in bin by bin\n", "arguments.update(dict(normBkgnd=1.0, fix_normBkgnd=True))\n", "# Remember to change \"twopar\" for two parameter fits\n", "arguments.update(dict(eps_eR_ee=0.0, fix_eps_eR_ee=True))\n", "#arguments.update(dict(eps_eR_ee=-0.472, error_eps_eR_ee=0.1))\n", "arguments.update(dict(eps_eL_ee=0.0, fix_eps_eL_ee=True))\n", "#arguments.update(dict(eps_eL_ee=-1.5, error_eps_eL_ee=0.1))\n", "\n", "arguments.update(dict(eps_eR_te=0.0, fix_eps_eR_te=True))\n", "arguments.update(dict(eps_eL_te=0.0, fix_eps_eL_te=True))\n", "arguments.update(dict(evis_res=0.0, fix_evis_res=True))\n", "arguments.update(dict(gL2=2.1785676, fix_gL2=True))\n", "arguments.update(dict(gR2=0.226576, fix_gR2=True))\n", "arguments.update(dict(gL=0.738, error_gL=0.1))\n", "arguments.update(dict(gR=0.238, error_gR=0.1))\n", "\n", "# Initialize a Minuit object. \n", "minuit = Minuit(fcn, **arguments)\n", "# Minimize fcn\n", "# You can put ncalls=10000 in the argument to migrad, but it's the default. \n", "migrad_result = minuit.migrad()\n", "# Display results. If you run this in Jupyter, it will have fancy HTML formatting.\n", "# This is not needed here, since we set print_level = 1 above. \n", "# minuit.print_fmin()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "name,value,err_parab,err_plus,err_minus\n", "gL 0.7381 0.2150 0.0533-0.2150\n", "normIBD 1.0005 0.1026 0.1031-0.1027\n", "norm 1.0000 0.0106 0.0106-0.0106\n", "gR 0.2383 0.6563 0.5476-0.6563\n" ] } ], "source": [ "# output Minuit variables and errors\n", "# Note that iminuit does not return a global correlation coefficient.\n", "unit11 = open(\"sinsqthw_nuebar_e_unit11.txt\",\"w\")\n", "line = 'name,value,err_parab,err_plus,err_minus'\n", "unit11.write(line + \"\\n\")\n", "print line\n", "\n", "# Run MINOS for every single variable. This returns a dictionary mapping each \n", "# variable name with the results of the MINOS fit.\n", "minos_results = minuit.minos()\n", "\n", "# Fetch migrad results into variables (actually, dictionaries).\n", "values = minuit.values\n", "errors = minuit.errors\n", "\n", "# Since we set print_level=1 above, the results of the MINOS fits will\n", "# be displayed. We also want to display the MINOS results in a file. \n", "\n", "for name in minos_results:\n", " minos_struct = minos_results[name]\n", " value = values[name]\n", " err_parab = errors[name]\n", " err_plus = minos_struct.upper\n", " err_minus = minos_struct.lower\n", " line = \"{:<11s}{:7.4f}{:7.4f}{:7.4f}{:7.4f}\".format(name,value,err_parab,err_plus,err_minus)\n", " unit11.write(line + \"\\n\")\n", " print line\n", " \n", "unit11.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n" ] } ], "source": [ "# option to make a contour\n", "if True:\n", " # contour variables\n", " npt=75\n", " xpt = np.zeros((npt))\n", " ypt = np.zeros((npt))\n", " sigmas = [1.,2.,3.]\n", " nsigs = len(sigmas)\n", " xcont = np.zeros((nsigs,npt+1))\n", " ycont = np.zeros((nsigs,npt+1))\n", " \n", " nsigs = 1\n", " for isigma in range(0,nsigs):\n", " nsigma = sigmas[isigma]\n", " print \"{:3.1f}\".format(nsigma)\n", " \n", " # Testing for equality with floating-point numbers is tricky, \n", " # so I'm checking if the difference between two numbers is small.\n", " errdef = 2.296 # default, for nsigma=1; sqrt(4.61) for 90% CL\n", " if abs(nsigma - 2.0) < 0.01:\n", " errdef = 6.18\n", " elif abs(nsigma - 3.0) < 0.01:\n", " errdef = 11.829\n", " elif abs(nsigma - 4.0) < 0.01:\n", " errdef = 19.334\n", " elif abs(nsigma - 5.0) < 0.01:\n", " errdef = 28.744 \n", " print 'SET ERRdef for Contour',errdef\n", " arguments['errordef'] = errdef\n", " \n", " # The contour is determined by letting each variable vary\n", " # by sigma*errdef**2\n", " sigma = 1.0\n", " mncontour_results = minuit.mncontour('gR','gL',npt,sigma)\n", " # mncontour_results = minuit.mncontour('eps_eR_ee','eps_eL_ee',npt,sigma)\n", " \n", " # mncontour returns a nested list of the form [[x1,y1]...[xn,yn]]\n", " nfound = len(mncontour_results)\n", " print 'npt,nfound: ',npt,nfound\n", " for i in range(0,nfound):\n", " xcont[isigma,i] = mncontour_results[i][0]\n", " ycont[isigma,i] = mncontour_results[i][1]\n", " xcont[isigma,nfound] = xcont[isigma,0]\n", " ycont[isigma,nfound] = ycont[isigma,0]\n", " \n", " for i in range(0,nfound+1):\n", " print \"{:3i}\".format(i),\n", " for j in range(0,nsigs):\n", " print \"{:9.5f}{:9.3f}\".format(xcont[j,i],ycont[j,i]),\n", " print \"\" # end of line" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }