{ "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)\n", " # eres = 0.02/sqrt(enu)\n", " #eres = 0.03/sqrt(enu) # <== This was used 3/13/18\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", " elen = sqrt((0.12/sqrt(enu))**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 = 15.000 MeV\n", "\n", "mass: 897.0 radius: 6.50 shield: 3.0 distance: 16.00\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.0640\n", "norm error = 0.050\n", "Eend error = 1e+06\n", " 1 0.01 0 2.343 1 0.002186 15\n", " 1 0.01122 0 1.861 1 0.002186 15\n", " 1 0.01259 0 1.479 1 0.002186 15\n", " 1 0.01413 0 1.175 1 0.002187 15\n", " 1 0.01585 0 0.9343 1 0.002187 15\n", " 1 0.01778 0 0.7427 1 0.002188 15\n", " 0.9684 0.01995 0 0.5905 1 0.002188 15\n", " 0.7701 0.02239 0 0.4696 1 0.002189 15\n", " 0.6126 0.02512 0 0.3736 1 0.002191 15\n", " 0.4875 0.02818 0 0.2973 1 0.002192 15\n", " 0.3882 0.03162 0 0.2367 1 0.002194 15\n", " 0.3092 0.03548 0 0.1886 1 0.002196 15\n", " 0.2466 0.03981 0 0.1503 1 0.002199 15\n", " 0.1968 0.04467 0 0.12 1 0.002203 15\n", " 0.1572 0.05012 0 0.09586 1 0.002208 15\n", " 0.1258 0.05623 0 0.07671 1 0.002214 15\n", " 0.1009 0.0631 0 0.0615 1 0.002222 15\n", " 0.08105 0.07079 0 0.04942 1 0.002231 15\n", " 0.06533 0.07943 0 0.03984 1 0.002244 15\n", " 0.05286 0.08913 0 0.03223 1 0.002259 15\n", " 0.04297 0.1 0 0.0262 1 0.00228 15\n", " 0.03513 0.1122 0 0.02142 1 0.002305 15\n", " 0.02893 0.1259 0 0.01764 1 0.002339 15\n", " 0.02404 0.1413 0 0.01466 1 0.002382 15\n", " 0.02021 0.1585 0 0.01232 1 0.002437 15\n", " 0.01721 0.1778 0 0.0105 1 0.002511 15\n", " 0.01491 0.1995 0 0.00909 1 0.002607 15\n", " 0.01317 0.2239 0 0.00803 1 0.002734 15\n", " 0.01191 0.2512 0 0.007261 1 0.002903 15\n", " 0.01106 0.2818 0 0.006741 1 0.003127 15\n", " 0.01057 0.3162 0 0.006444 1 0.003423 15\n", " 0.01041 0.3548 0 0.006349 1 0.003807 15\n", " 0.01056 0.3981 0 0.006438 1 0.004289 15\n", " 0.01094 0.4467 0 0.006669 1 0.00485 15\n", " 0.01133 0.5012 0 0.006907 1 0.005375 15\n", " 0.01125 0.5623 0 0.00686 1 0.005587 15\n", " 0.01036 0.631 0 0.006316 1 0.005249 15\n", " 0.009029 0.7079 0 0.005506 1 0.00454 15\n", " 0.007858 0.7943 0 0.004792 1 0.003801 15\n", " 0.007115 0.8913 0 0.004338 1 0.003208 15\n", " 0.006815 1 0 0.004155 1 0.002786 15\n", " 0.00687 1.122 0 0.004189 1 0.002509 15\n", " 0.007115 1.259 0 0.004339 1 0.00235 15\n", " 0.007285 1.413 0 0.004442 1 0.002288 15\n", " 0.007228 1.585 0 0.004407 1 0.00231 15\n", " 0.007093 1.778 0 0.004325 1 0.002401 15\n", " 0.007075 1.995 0 0.004314 1 0.002535 15\n", " 0.007206 2.239 0 0.004394 1 0.002663 15\n", " 0.007398 2.512 0 0.004511 1 0.002741 15\n", " 0.007595 2.818 0 0.004631 1 0.002774 15\n", " 0.007822 3.162 0 0.00477 1 0.002803 15\n", " 0.008118 3.548 0 0.00495 1 0.002854 15\n", " 0.008513 3.981 0 0.005191 1 0.002941 15\n", " 0.009015 4.467 0 0.005497 1 0.003081 15\n", " 0.009656 5.012 0 0.005888 1 0.003271 15\n", " 0.0105 5.623 0 0.006401 1 0.00351 15\n", " 0.01162 6.31 0 0.007086 1 0.00382 15\n", " 0.01314 7.079 0 0.008015 1 0.004248 15\n", " 0.01527 7.943 0 0.009313 1 0.004866 15\n", " 0.01834 8.913 0 0.01118 1 0.005768 15\n", " 0.02295 10 0 0.014 1 0.007137 15\n", " 0.03027 11.22 0 0.01846 1 0.009331 15\n", " 0.04267 12.59 0 0.02602 1 0.01308 15\n", " 0.0653 14.13 0 0.03982 1 0.01995 15\n", " 0.1063 15.85 0 0.06481 1 0.03241 15\n", " 0.1533 17.78 0 0.09349 1 0.04673 15\n", " 0.161 19.95 0 0.0982 1 0.04908 15\n", " 0.154 22.39 0 0.09388 1 0.04693 15\n", " 0.1562 25.12 0 0.09526 1 0.04762 15\n", " 0.1671 28.18 0 0.1019 1 0.05094 15\n", " 0.1794 31.62 0 0.1094 1 0.05467 15\n", " 0.185 35.48 0 0.1128 1 0.05638 15\n", " 0.1844 39.81 0 0.1124 1 0.05619 15\n", " 0.1837 44.67 0 0.112 1 0.05597 15\n", " 0.1852 50.12 0 0.1129 1 0.05644 15\n", " 0.1868 56.23 0 0.1139 1 0.05692 15\n", " 0.1871 63.1 0 0.1141 1 0.057 15\n", " 0.1871 70.79 0 0.1141 1 0.05703 15\n", " 0.1872 79.43 0 0.1142 1 0.05706 15\n", " 0.1873 89.13 0 0.1142 1 0.05707 15\n", " 0.1873 100 0 0.1142 1 0.05708 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.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.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() # change -1 to +1\n", " # enu,levt,n_nuebar_IBD,n_nuebar\n", " fields = line.split()\n", "\n", "# for i in range(0,4):\n", "# print float(fields[i])\n", "\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 }