{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "This is the python equivalent of IBD_disapp.f written by Mike Shaevitz. It was adapted for python by William Seligman. \n", "\n", "_Markdown cells (like this one) indicate a note added by Seligman during his conversion process._\n", "\n", "_First note: The odds that I've made an arithmetic error in the code conversion is probably 100%._\n", "\n", "__Second note: This program requires iminuit. This is not installed by default in most Python distributions. The following command might do this (check with your sysadmin):__\n", "\n", "``pip install iminuit``\n", "\n", "The iminuit documentation is [here](http://iminuit.readthedocs.io/en/latest/api.html#api-doc).\n", "\n", "_It's standard in Jupyter to develop code in cells: Make sure the code in a cell is working, then move onto execute the next one. To execute all the cells as if they were a single program, select **Cell->Run All** in the menu above._" ] }, { "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", "# \"isnan(x)\" tests whether 'x' is NaN (not a number); e.g., from dividing by zero. \n", "from math import sqrt,isnan\n", "from random import random" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# For the sine integral\n", "import scipy\n", "from scipy.special import sici\n", "# For sine and cosine. (These are available without the following statement,\n", "# but you'd have to remember to use math.sin() and math.cos().)\n", "from math import sin, cos" ] }, { "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": 3, "metadata": {}, "outputs": [], "source": [ "# L/E bins and limits\n", "global nle\n", "global ile\n", "global nev\n", "global dle_max\n", "global dle\n", "global dat_array\n", "global levt_smr\n", "global evt_le" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the osc_int function" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def osc_int(dm2,e1,e2,l1,l2,esigma):\n", "\n", " # this routine calculates the average of the \n", " # oscillation formula (sin(1.27*dm2*l/e)**2)\n", " # over the distance range l1 to l2 and energy range e1 to e2\n", " # (The code for averages separately over e and l are included as comments)\n", " \n", " # formula to average over both e1 to e2 and l1 to l2\n", " k = 1.27*dm2\n", " a11 = 2.*k*l1/e1\n", " a12 = 2.*k*l1/e2\n", " a21 = 2.*k*l2/e1\n", " a22 = 2.*k*l2/e2\n", " # Python does not have the CERNLIB DSININ function (sine integral).\n", " # Instead, we must use the SciPy function sici which returns both the\n", " # sine and cosine integrals as a two-element list, then pick the one we want to use. \n", " (si_a11,ci_a11) = sici(a11)\n", " (si_a12,ci_a12) = sici(a12)\n", " (si_a21,ci_a21) = sici(a21)\n", " (si_a22,ci_a22) = sici(a22)\n", " posc_full = -4.*e2*k*l2-4*k*k*si_a21*l2*l2 \\\n", " +4.*k*k*si_a22*l2*l2 \\\n", " +4.*e1*k*l2-sin(a21)*e1*e1+sin(a22)*e2*e2 \\\n", " -2.*k*cos(a21)*e1*l2+2.*k*cos(a22)*e2*l2 \\\n", " +4.*e2*k*l1+4.*k*k*si_a11*l1*l1 \\\n", " -4.*k*k*si_a12*l1*l1-4.*e1*k*l1+sin(a11)*e1*e1 \\\n", " -sin(a12)*e2*e2+2.*k*cos(a11)*e1*l1 \\\n", " -2.*k*cos(a12)*e2*l1\n", " posc_full = -posc_full/(8.*k)/(e2-e1)/(l2-l1)\n", " \n", " # probability should be 0.5 if oscillations less than energy resolution\n", " #eavg = (e1+e2)/2.\n", " #lavg = (l1+l2)/2.\n", " #n = 1.27*dm2*lavg/(eavg*3.14159)\n", " #enext = 1.27*dm2*lavg/((n+0.5)*3.14159)\n", " # if (n > 2 and abs(eavg-enext)/eavg < esigma:\n", " # posc_full = 0.5\n", " \n", " # isnan means \"is NaN\"; in other words, the calculation has a divide by zero somewhere.\n", " if isnan(posc_full):\n", " print \"Debug! posc NaN k=\",k,\" si_a11=\",si_a11,\" si_a12=\",si_a12, \\\n", " \"si_a21=\",si_a21,\" si_a22=\",si_a22\n", " print \"Debug! a11=\",a11,\" a12=\",a12,\" a21=\",a21,\" a22=\",a22\n", " print \"Debug! dm2=\",dm2,\" e2=\",e2,\" e1=\",e1,\" l2=\",l2,\" l1=\",l1\n", " \n", " # output value\n", " return posc_full" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Define the binLE function for use by the main program._" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def binLE(dm2,sinsq2th,Eend,evt_le_output,fullosc_le_output):\n", " \n", " # In the FORTRAN code, these variables were passed as arguments. \n", " # For Python, we use the global statement to access their definitions\n", " # in the main program. \n", " # L/E bins and limits\n", " global dat_array\n", " global levt_smr\n", " global nle\n", " global dle_max\n", " global nev\n", " \n", " # Kopp_Maltoni_Schwetz 1103.4570 - Ue4 and Ue5\n", " # dm2_41 = 0.47\n", " # U_e4 = 0.384\n", " # dm2_51 = 0.87\n", " # U_e5 = 0.414\n", " dm2_41 = 0.47\n", " U_e4 = 0.128\n", " dm2_51 = 0.87\n", " U_e5 = 0.138\n", " \n", " # bin events and fullosc in L/E for this dm2\n", " dle = dle_max/nle\n", " for iev in range(0,nev):\n", " enu = dat_array[0,iev]\n", " levt = dat_array[1,iev]\n", " spec_corr = \\\n", " ((enu*(Eend-enu))/(enu*(15.-enu)))**2 \\\n", " /((8.*(Eend-8.))/(8.*(15.-8.)))**2\n", " n_nuebar_IBD = dat_array[2,iev]*spec_corr\n", " ile = int(levt/enu/dle)\n", " ile = min(nle-1,max(0,ile))\n", " \n", " # 6.4%/Sqrt(E) and ~12 cm /Sqrt(E) - Lindley email 3/1/12\n", " # Use x1.73 to go from sigma to Half_width@base\n", " # interaction position has flat dist from 0 to 175cm\n", " # with 11B blanket, production position is 100 cm\n", " # eres = 0.064/sqrt(enu) # nominal for KamLAND\n", " # eres = 0.02/sqrt(enu)\n", " #eres = 0.03/sqrt(enu) # <== This was used 3/13/18\n", " eres = (-0.0839+0.349*sqrt(max(enu-1.8,0.5))+0.0397*max(enu-1.8,0.5))/enu # SuperK (aXiv:1606.07538)\n", " \n", " # Python doesn't have static initialization values the way\n", " # FORTRAN does. The following tests if the variable 'firstLE'\n", " # exists. If it doesn't, after the print statement the\n", " # variable is defined globally so it will exist for subsequent\n", " # executions for this routine. \n", " global firstLE\n", " try:\n", " firstLE\n", " except NameError:\n", " print (\"energy resolution = {:7.4f}\".format(eres*sqrt(enu)))\n", " firstLE = 1\n", " #vert_res = 0.12 # nominal KamLAND in m\n", " vert_res = 1.70 # SuperK approx vertex resoluion (aXiv:1606.07538)\n", " elen = sqrt((vert_res/sqrt(max(enu-1.8,0.5)))**2+(0.41)**2)\n", " l = levt\n", " # use box smearing\n", " # elen = 0.001\n", " # l = levt_smr[iev]\n", " # oscprob 3+1\n", " oscprob31 = True\n", " if oscprob31:\n", " oscprob = osc_int(dm2,enu*(1.-eres*1.73),\n", " enu*(1.+eres*1.73),l-elen*1.73,l+elen*1.73,eres)\n", " if isnan(oscprob):\n", " print \"Debug! dm2=\",dm2,\" enu=\",enu,\" eres=\",eres,\" elen=\",elen\n", " \n", " # oscprob 3+2 (Kopp/Maltoni/Schwetz)\n", " oscprob32 = False\n", " if oscprob32:\n", " osc41 = osc_int(dm2_41,enu*(1.-eres*1.73),\n", " enu*(1.+eres*1.73),levt-elen*1.73,levt+elen*1.73,eres)\n", " osc51 = osc_int(dm2_51,enu*(1.-eres*1.73),\n", " enu*(1.+eres*1.73),levt-elen*1.73,levt+elen*1.73,eres)\n", " osc54 = osc_int(dm2_51-dm2_41,enu*(1.-eres*1.73),\n", " enu*(1.+eres*1.73),levt-elen*1.73,levt+elen*1.73,eres)\n", " oscprob = 4.*((1.-U_e4**2-U_e5**2)*\n", " (U_e4**2*osc41+U_e5**2*osc51)+(U_e4**2*U_e5**2*osc54))\n", " \n", " if levt/enu <= dle_max:\n", " evt_le_output[ile] += n_nuebar_IBD*(1.-oscprob*sinsq2th)\n", " fullosc_le_output[ile] += n_nuebar_IBD*oscprob\n", " \n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the function to be minimized by Minuit. Note that the function is defined in terms of _f(x,y,z,...)_ instead of the FORTRAN-based method of passing parameter arrays. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def fcn(sinsq2th,norm,dm2,Eend):\n", " \n", " # L/E bins and limits\n", " # from the old FORTRAN COMMON block\n", " # integer nle,ile\n", " # real*8 dle_max,dle,enorm\n", " # real*8 evt_le(10000),fullosc_le(10000)\n", " # common /lebins/ dle_max,dle,evt_le,fullosc_le,nle,ile\n", " global evt_le\n", " \n", " # pull term contribution\n", " enorm = 0.05\n", " chi2 = (1.-norm)**2/(enorm)**2\n", " eEend = 1.e6 # MeV\n", " chi2 = chi2 + (Eend-15.0)**2/(eEend)**2\n", " # Python doesn't have static initialization values the way\n", " # FORTRAN does. The following tests if the variable 'firstFCN'\n", " # exists. If it doesn't, after the print statement the\n", " # variable is defined globally so it will exist for subsequent\n", " # executions for this routine. \n", " global firstFCN\n", " try:\n", " firstFCN\n", " except NameError:\n", " print (\"norm error = {:7.3f}\".format(enorm))\n", " print (\"Eend error = {:12.4g}\".format(eEend))\n", " firstFCN = 1\n", " \n", " # bin events in L/E bins\n", " evt_le_pred = np.zeros((nle))\n", " fullosc_le_pred = np.zeros((nle))\n", " binLE(dm2,sinsq2th,Eend,evt_le_pred,fullosc_le_pred)\n", "\n", " # calculate chi2 over L/E bins\n", " for ile in range (0,nle):\n", " \n", " # pred = norm*(evt_le[ile]-sinsq2th*fullosc_le[ile])\n", " \n", " pred = norm*evt_le_pred[ile]\n", " \n", " data = evt_le[ile]\n", " \n", " if isnan(pred) or isnan(data):\n", " print \"Zingo! ile=\",ile,\" pred=\",pred,\" data=\",data, \\\n", " \" chi2=\",chi2,\" evt_le_pred[ile]=\",evt_le_pred[ile],\" norm=\",norm\n", "\n", " # if(data > 1.): # use with likelihood statistic\n", " if data > 10.: # use with chi2 statistic\n", " # use chi2 statistic\n", " chi2 = chi2 + \\\n", " (data-pred)**2 \\\n", " /max(1.,data)\n", " # use binned max likelihood statistic\n", " # 2*(pred-data-data*log(pred/data))\n", "\n", " # if chi2 == 0.:\n", " # print \"Zingo! chi2=0, sinsq2th=\",sinsq2th,\" norm=\",norm,\" dm2=\",dm2,\" Eend=\",Eend\n", " \n", " return chi2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main routine" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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 (MeV): 15.000 eff_det: 0.690 eff_Gd_tag: 0.690\n", "\n", "mass: 22500.0 height: 32.20 radius: 14.91 shield: 9.25 distance: 9.25\n", "\n", "nyrs_dat_array: 4.00 nyrs run: 2.40\n", "elimit_factor = 1.640\n", "sinsq2th_limit dm2 sinsq2th esinsq2th fnorm efnorm Eend eEend\n", "Smearing box size: +/- 1.00 m\n", "energy resolution = 0.6672\n", "norm error = 0.050\n", "Eend error = 1e+06\n", " 0.6473 0.01 0 0.3947 1 0.0009409 15\n", " 0.5146 0.01122 0 0.3138 1 0.0009411 15\n", " 0.4092 0.01259 0 0.2495 1 0.0009414 15\n", " 0.3255 0.01413 0 0.1985 1 0.0009417 15\n", " 0.2589 0.01585 0 0.1579 1 0.0009421 15\n", " 0.2061 0.01778 0 0.1257 1 0.0009426 15\n", " 0.1641 0.01995 0 0.1001 1 0.0009433 15\n", " 0.1308 0.02239 0 0.07976 1 0.0009441 15\n", " 0.1043 0.02512 0 0.06361 1 0.0009452 15\n", " 0.08329 0.02818 0 0.05079 1 0.0009465 15\n", " 0.06659 0.03162 0 0.0406 1 0.0009482 15\n", " 0.05332 0.03548 0 0.03251 1 0.0009503 15\n", " 0.04279 0.03981 0 0.02609 1 0.000953 15\n", " 0.03443 0.04467 0 0.02099 1 0.0009564 15\n", " 0.02779 0.05012 0 0.01694 1 0.0009608 15\n", " 0.02252 0.05623 0 0.01373 1 0.0009663 15\n", " 0.01834 0.0631 0 0.01118 1 0.0009734 15\n", " 0.01503 0.07079 0 0.009167 1 0.0009824 15\n", " 0.01242 0.07943 0 0.007572 1 0.000994 15\n", " 0.01036 0.08913 0 0.006315 1 0.001009 15\n", " 0.008736 0.1 0 0.005327 1 0.001028 15\n", " 0.007474 0.1122 0 0.004557 1 0.001054 15\n", " 0.0065 0.1259 0 0.003963 1 0.001086 15\n", " 0.005761 0.1413 0 0.003513 1 0.001129 15\n", " 0.005218 0.1585 0 0.003182 1 0.001184 15\n", " 0.004837 0.1778 0 0.00295 1 0.001256 15\n", " 0.004595 0.1995 0 0.002802 1 0.001346 15\n", " 0.004474 0.2239 0 0.002728 1 0.001459 15\n", " 0.004467 0.2512 0 0.002724 1 0.001598 15\n", " 0.00457 0.2818 0 0.002786 1 0.001766 15\n", " 0.004764 0.3162 0 0.002905 1 0.00196 15\n", " 0.004986 0.3548 0 0.00304 1 0.002148 15\n", " 0.005111 0.3981 0 0.003117 1 0.00227 15\n", " 0.005026 0.4467 0 0.003064 1 0.002264 15\n", " 0.004756 0.5012 0 0.0029 1 0.002139 15\n", " 0.004454 0.5623 0 0.002716 1 0.001969 15\n", " 0.004244 0.631 0 0.002588 1 0.001817 15\n", " 0.004181 0.7079 0 0.002549 1 0.00171 15\n", " 0.004272 0.7943 0 0.002605 1 0.001649 15\n", " 0.004515 0.8913 0 0.002753 1 0.001631 15\n", " 0.004905 1 0 0.002991 1 0.001656 15\n", " 0.005433 1.122 0 0.003313 1 0.001727 15\n", " 0.006057 1.259 0 0.003693 1 0.001849 15\n", " 0.006722 1.413 0 0.004099 1 0.002021 15\n", " 0.007452 1.585 0 0.004544 1 0.002254 15\n", " 0.008363 1.778 0 0.005099 1 0.00257 15\n", " 0.009581 1.995 0 0.005842 1 0.002984 15\n", " 0.01124 2.239 0 0.006856 1 0.00352 15\n", " 0.01347 2.512 0 0.008213 1 0.004206 15\n", " 0.01629 2.818 0 0.009931 1 0.00505 15\n", " 0.01952 3.162 0 0.0119 1 0.00601 15\n", " 0.02311 3.548 0 0.01409 1 0.007081 15\n", " 0.02741 3.981 0 0.01672 1 0.008379 15\n", " 0.03271 4.467 0 0.01995 1 0.009987 15\n", " 0.03928 5.012 0 0.02395 1 0.01199 15\n", " 0.04813 5.623 0 0.02935 1 0.01469 15\n", " 0.06141 6.31 0 0.03744 1 0.01874 15\n", " 0.08169 7.079 0 0.04981 1 0.02492 15\n", " 0.1136 7.943 0 0.06928 1 0.03464 15\n", " 0.1684 8.913 0 0.1027 1 0.05133 15\n", " 0.2692 10 0 0.1641 1 0.08206 15\n", " 0.4396 11.22 0 0.268 1 0.134 15\n", " 0.6222 12.59 0 0.3794 1 0.1897 15\n", " 0.6303 14.13 0 0.3844 1 0.1922 15\n", " 0.6075 15.85 0 0.3704 1 0.1852 15\n", " 0.645 17.78 0 0.3933 1 0.1966 15\n", " 0.7581 19.95 0 0.4623 1 0.2311 15\n", " 0.9572 22.39 0 0.5837 1 0.2918 15\n", " 0.03605 25.12 0 0.02198 1 0.01099 15\n", " 0.03605 28.18 0 0.02198 1 0.01099 15\n", " 0.03605 31.62 0 0.02198 1 0.01099 15\n", " 0.03605 35.48 0 0.02198 1 0.01099 15\n", " 0.03605 39.81 0 0.02198 1 0.01099 15\n", " 0.03605 44.67 0 0.02198 1 0.01099 15\n", " 0.03605 50.12 0 0.02198 1 0.01099 15\n", " 0.03605 56.23 0 0.02198 1 0.01099 15\n", " 0.03605 63.1 0 0.02198 1 0.01099 15\n", " 0.03605 70.79 0 0.02198 1 0.01099 15\n", " 0.03605 79.43 0 0.02198 1 0.01099 15\n", " 0.03605 89.13 0 0.02198 1 0.01099 15\n", " 0.03605 100 0 0.02198 1 0.01099 15\n" ] } ], "source": [ "# The contents of the old FORTRAN COMMON block.\n", "# Note that these are shared via global statements in a previous cell.\n", "# Also note that we can determine array or list sizes dynamically,\n", "# so we don't need to set artificial limits like 200000000.\n", "\n", "# L/E bins and limits\n", "# integer nle,ile,nev\n", "# real*8 dle_max,dle\n", "# real*8 evt_le(10000),fullosc_le(10000)\n", "# real*8 dat_array(4,200000000)\n", "# real*8 levt_smr(200000000)\n", "# common /lebins/ nev,dle_max,dle,nle,ile,\n", "# > evt_le,fullosc_le,dat_array,levt_smr\n", "\n", "# Kopp_Maltoni_Schwetz 1103.4570 - Ue4 and Ue5\n", "# dm2_41 = 0.47\n", "# U_e4 = 0.384\n", "# dm2_51 = 0.87\n", "# U_e5 = 0.414\n", "dm2_41 = 0.47\n", "U_e4 = 0.128\n", "dm2_51 = 0.87\n", "U_e5 = 0.138\n", "\n", "# setup\n", "# nyrs = 3.0*0.8/1.68 #<== 1.68 for 12.5m equal stats\n", "nyrs = 3.0*0.8\n", "elimit_factor = 1.64 # 95% CL one-sided\n", "# nle = 400 # baseline\n", "nle = 1000 # expanded - standard sensitivity run\n", "# nle = 60 # <= standard L/E plot\n", "# nle = 200\n", "# nle = 30\n", "dle_max = 10.\n", "\n", "# print out run parameters\n", "with open(\"current_istope_event_calc_run_SuperK.txt\",\"r\") as f:\n", " contents = f.readlines() # Read the file into a list of lines.\n", " \n", "for i in range(0,5):\n", " print contents[i]\n", "\n", "# The numbers we want are in the second line of the file. \n", "line = contents[1].strip()\n", "fields = line.split()\n", "# This comes from looking at the line printed out above:\n", "# \"total flux: 1.147e+23 nyrs(with DF): 4.00 protons/yr: 1.97e+24\"\n", "# We want the 6th non-blank field (I think).\n", "nyrs_dat_array = float(fields[5]) \n", "\n", "print \"nyrs_dat_array: {:5.2f} nyrs run: {:5.2f}\".format(nyrs_dat_array,nyrs)\n", "print \"elimit_factor = {:6.3f}\".format(elimit_factor)\n", "print 'sinsq2th_limit dm2 sinsq2th ' \\\n", " 'esinsq2th fnorm efnorm Eend eEend'\n", "\n", "# read in data array\n", "with open(\"nuebar_event_array_SuperK.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", "# This comes from looking at the line in the file:\n", "# \"10000 0.238\"\n", "# We want the 1st non-blank field.\n", "nev = int(fields[0])\n", "\n", "# box region for smearing distance\n", "box = 1.0 # m - need to remove comment on levt_smr in binLE also\n", "\n", "levt_smr = np.zeros((nev))\n", "dat_array = np.zeros((4,nev)) # Define 2D array\n", "for iev in range(0,nev):\n", " line = contents[iev+1].strip()\n", " # enu,levt,n_nuebar_IBD,n_nuebar\n", " fields = line.split()\n", " for i in range(0,4):\n", " dat_array[i,iev] = float(fields[i])\n", " # correct for years\n", " dat_array[2,iev] *= nyrs/nyrs_dat_array\n", " # throw production point in 1m box\n", " # \"random()\" with no arguments generates a random number between 0 and 1. \n", " x = (1.-2.*random())*box\n", " y = (1.-2.*random())*box\n", " levt = dat_array[1,iev]\n", " z = levt + (1.-2.*random())*box\n", " levt_smr[iev] = sqrt(x**2+y**2+z**2)\n", "\n", "print \"Smearing box size: +/- {:6.2f} m\".format(box)\n", "\n", "# dm2 loop\n", "for idm2 in range(0,81):\n", "# for idm2 in range(40,41):\n", " dm2 = 10**(-2.+float(idm2)/20.)\n", " \n", " # bin events in L/E bins\n", " sinsq2th0 = 0.0\n", " evt_le = np.zeros((nle))\n", " fullosc_le = np.zeros((nle))\n", " binLE(dm2,sinsq2th0,15.0,evt_le,fullosc_le)\n", " \n", " # initialize minuit\n", " # ===================================================== \n", " # In Python, we fiddle with the operation of Minuit by\n", " # passing keyword arguments to the Minuit method. Because there are\n", " # so many arguments for this code, we're using a Python\n", " # dictionary to set the values. \n", " # =====================================================\n", " # set up parameters that vary\n", " arguments = dict(sinsq2th=0.0, error_sinsq2th=0.01)\n", " arguments.update(dict(norm=1.0, error_norm=0.1))\n", " arguments.update(dict(dm2=dm2, fix_dm2=True))\n", " arguments.update(dict(Eend=15.0, fix_Eend=True))\n", " # set the print level\n", " # Set to 1 for output at end of migrad/hesse/minos \n", " arguments.update(dict(print_level=0))\n", " # set value for error calculation\n", " arguments.update(dict(errordef=1.0))\n", "\n", " # Initialize a Minuit object. \n", " minuit = Minuit(fcn, **arguments)\n", " # Minimize fcn(sinsq2th,norm,dm2,Eend)\n", " # You can put ncalls=10000 in the argument to migrad, but it's the default. \n", " result = minuit.migrad()\n", " \n", " # Display results. If you run this in Jupyter, it will have fancy HTML formatting. \n", " # minuit.print_fmin() # <== This prints out all the results of the fit - very nicely.\n", " # minuit.print_matrix()\n", " \n", " # If you wanted to run minos for parameter sinsq2th (for example), it would be:\n", " # result.minos(sinsq2th)\n", "\n", " # Fetch results into variables (actually, dictionaries).\n", " values = minuit.values\n", " errors = minuit.errors\n", " sinsq2th = values['sinsq2th']\n", " esinsq2th = errors['sinsq2th']\n", " fnorm = values['norm']\n", " efnorm = errors['norm']\n", " Eend = values['Eend']\n", " eEend = errors['Eend']\n", " \n", " formats = \"\" # construct the equivalent of FORTRAN's 8G12.4\n", " for i in range(1,8):\n", " formats += \"{:12.4g}\"\n", " print formats.format(min(1.,esinsq2th*elimit_factor), \\\n", " dm2,sinsq2th,esinsq2th,fnorm,efnorm,Eend,eEend)\n", " \n", " # ==========================================================\n", " if False:\n", " # sinsq2th_plot = min(1.,esinsq2th*1.28)\n", " sinsq2th_plot = 0.1\n", " print 'L/E Result for dm2,sinsq2th: {:8.4f}{:8.4f}' \\\n", " .format(dm2,sinsq2th_plot)\n", " print 'ile,le,evt_le,fullosc_le,prob_evts,' \\\n", " ,'eprob,oscprob,oscprob_bin'\n", " for ile in range (0,nle-1):\n", " # oscprob for 3+1\n", " if True:\n", " oscprob_bin = sinsq2th_plot*sin(1.27*dm2*(ile+0.5)*dle)**2\n", " evt_meas = evt_le[ile]-sinsq2th_plot*fullosc_le[ile]\n", " # oscprob for 3+2\n", " if False:\n", " osc41 = sin(1.27*dm2_41*(ile-0.5)*dle)**2\n", " osc51 = sin(1.27*dm2_51*(ile-0.5)*dle)**2\n", " osc54 = sin(1.27*(dm2_51-dm2_41)*(ile+0.5)*dle)**2\n", " oscprob_bin = 4.*((1.-U_e4**2-U_e5**2)* \\\n", " (U_e4**2*osc41+U_e5**2*osc51)+(U_e4**2*U_e5**2*osc54))\n", " evt_meas = evt_le[ile]-fullosc_le[ile]\n", " \n", " formats = \"{:3d}\" # construct the equivalent of FORTRAN's I3,8G12.4\n", " for i in range(1,8):\n", " formats += \"{:12.4g}\"\n", " print formats.format(ile,(ile+0.5)*dle,evt_le[ile],fullosc_le[ile] \\\n", " ,evt_meas/max(1.,evt_le[ile]) \\\n", " ,1./sqrt(max(1.,evt_le[ile])) \\\n", " ,sinsq2th_plot*fullosc_le[ile]/max(1.,evt_le[ile]) \\\n", " ,oscprob_bin)" ] }, { "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 }