{ "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: 100.0 radius: 3.01 shield: 5.0 distance: 14.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 13.19 1 0.008268 15\n", " 1 0.01122 0 10.48 1 0.008268 15\n", " 1 0.01259 0 8.327 1 0.008269 15\n", " 1 0.01413 0 6.616 1 0.00827 15\n", " 1 0.01585 0 5.257 1 0.008271 15\n", " 1 0.01778 0 4.177 1 0.008272 15\n", " 1 0.01995 0 3.32 1 0.008273 15\n", " 1 0.02239 0 2.638 1 0.008275 15\n", " 1 0.02512 0 2.097 1 0.008278 15\n", " 1 0.02818 0 1.668 1 0.008281 15\n", " 1 0.03162 0 1.326 1 0.008285 15\n", " 1 0.03548 0 1.055 1 0.00829 15\n", " 1 0.03981 0 0.8396 1 0.008296 15\n", " 1 0.04467 0 0.6685 1 0.008304 15\n", " 0.8735 0.05012 0 0.5326 1 0.008314 15\n", " 0.6965 0.05623 0 0.4247 1 0.008326 15\n", " 0.5559 0.0631 0 0.3389 1 0.008342 15\n", " 0.4442 0.07079 0 0.2709 1 0.008362 15\n", " 0.3555 0.07943 0 0.2168 1 0.008388 15\n", " 0.2851 0.08913 0 0.1738 1 0.00842 15\n", " 0.2292 0.1 0 0.1398 1 0.008461 15\n", " 0.1848 0.1122 0 0.1127 1 0.008513 15\n", " 0.1497 0.1259 0 0.09125 1 0.00858 15\n", " 0.1218 0.1413 0 0.07425 1 0.008666 15\n", " 0.09972 0.1585 0 0.06081 1 0.008778 15\n", " 0.08232 0.1778 0 0.05019 1 0.008923 15\n", " 0.06865 0.1995 0 0.04186 1 0.009113 15\n", " 0.058 0.2239 0 0.03537 1 0.009366 15\n", " 0.04982 0.2512 0 0.03038 1 0.009706 15\n", " 0.0437 0.2818 0 0.02665 1 0.01017 15\n", " 0.03938 0.3162 0 0.02402 1 0.01082 15\n", " 0.03673 0.3548 0 0.0224 1 0.01174 15\n", " 0.03574 0.3981 0 0.02179 1 0.01308 15\n", " 0.03656 0.4467 0 0.02229 1 0.01504 15\n", " 0.03916 0.5012 0 0.02388 1 0.01773 15\n", " 0.04204 0.5623 0 0.02563 1 0.0205 15\n", " 0.04085 0.631 0 0.02491 1 0.02094 15\n", " 0.03477 0.7079 0 0.0212 1 0.0182 15\n", " 0.02866 0.7943 0 0.01748 1 0.01483 15\n", " 0.02465 0.8913 0 0.01503 1 0.01216 15\n", " 0.02257 1 0 0.01376 1 0.01019 15\n", " 0.02219 1.122 0 0.01353 1 0.008816 15\n", " 0.02331 1.259 0 0.01421 1 0.007878 15\n", " 0.02486 1.413 0 0.01516 1 0.007238 15\n", " 0.02486 1.585 0 0.01516 1 0.006893 15\n", " 0.02313 1.778 0 0.0141 1 0.006939 15\n", " 0.02223 1.995 0 0.01355 1 0.007502 15\n", " 0.02273 2.239 0 0.01386 1 0.00845 15\n", " 0.02362 2.512 0 0.0144 1 0.009213 15\n", " 0.02381 2.818 0 0.01452 1 0.009125 15\n", " 0.02418 3.162 0 0.01474 1 0.008718 15\n", " 0.02519 3.548 0 0.01536 1 0.008589 15\n", " 0.02628 3.981 0 0.01602 1 0.008924 15\n", " 0.02755 4.467 0 0.0168 1 0.009528 15\n", " 0.02933 5.012 0 0.01789 1 0.01008 15\n", " 0.03166 5.623 0 0.0193 1 0.0106 15\n", " 0.03466 6.31 0 0.02114 1 0.01138 15\n", " 0.03872 7.079 0 0.02361 1 0.01256 15\n", " 0.04414 7.943 0 0.02692 1 0.01412 15\n", " 0.05172 8.913 0 0.03154 1 0.01628 15\n", " 0.0624 10 0 0.03805 1 0.0194 15\n", " 0.07752 11.22 0 0.04727 1 0.02389 15\n", " 0.09848 12.59 0 0.06005 1 0.03012 15\n", " 0.1244 14.13 0 0.07587 1 0.03791 15\n", " 0.1485 15.85 0 0.09056 1 0.04514 15\n", " 0.1614 17.78 0 0.09843 1 0.04902 15\n", " 0.1641 19.95 0 0.1 1 0.04982 15\n", " 0.1636 22.39 0 0.09977 1 0.04969 15\n", " 0.164 25.12 0 0.09999 1 0.04979 15\n", " 0.1651 28.18 0 0.1007 1 0.05013 15\n", " 0.166 31.62 0 0.1012 1 0.0504 15\n", " 0.1663 35.48 0 0.1014 1 0.05049 15\n", " 0.1662 39.81 0 0.1014 1 0.05047 15\n", " 0.1661 44.67 0 0.1013 1 0.05045 15\n", " 0.1663 50.12 0 0.1014 1 0.05049 15\n", " 0.1665 56.23 0 0.1015 1 0.05055 15\n", " 0.1666 63.1 0 0.1016 1 0.05057 15\n", " 0.1666 70.79 0 0.1016 1 0.05058 15\n", " 0.1666 79.43 0 0.1016 1 0.05058 15\n", " 0.1666 89.13 0 0.1016 1 0.05058 15\n", " 0.1666 100 0 0.1016 1 0.05058 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_Borexino.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_Borexino.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 }