{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "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": 71, "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": 72, "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": 73, "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": 74, "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": 75, "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 dle\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.10/sqrt(enu) # <== CHANDLER\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", " elen = sqrt((0.1)**2+(0.41)**2) # <== CHANDLER??\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": 76, "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", " #print 'fcn call: ', chi2,sinsq2th,norm,dm2,Eend\n", " \n", " return chi2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main routine" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "24200 1210000 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: 40.5 shield: 2.8 distance source to detector: 4.00 levt_max 12.894\n", "\n", "nyrs_dat_array: 4.00 nyrs run: 4.00\n", "elimit_factor = 1.640\n", "sinsq2th_limit dm2 sinsq2th esinsq2th fnorm efnorm Eend eEend\n", "nev: 24200\n", "Smearing box size: +/- 1.00 m\n", " 0.008743 1 0 0.005331 1 0.003939 15\n", "L/E Result for dm2,sinsq2th: 1.0000 0.1000\n", "ile,le,evt_le,fullosc_le,prob_evts, eprob,oscprob,oscprob_bin\n", " 0 0.02 0 0 0 1 0 6.45e-05\n", " 1 0.06 0 0 0 1 0 0.0005795\n", " 2 0.1 0 0 0 1 0 0.001604\n", " 3 0.14 0 0 0 1 0 0.003128\n", " 4 0.18 0 0 0 1 0 0.005135\n", " 5 0.22 0 0 0 1 0 0.007605\n", " 6 0.26 0.1065 0.01282 0.1052 1 0.001282 0.01051\n", " 7 0.3 162.5 24.5 0.9849 0.07845 0.01508 0.01383\n", " 8 0.34 2088 374.3 0.9821 0.02188 0.01792 0.01751\n", " 9 0.38 8579 1899 0.9779 0.0108 0.02213 0.02154\n", " 10 0.42 1.273e+04 3334 0.9738 0.008864 0.0262 0.02585\n", " 11 0.46 1.593e+04 4897 0.9693 0.007924 0.03075 0.03042\n", " 12 0.5 2.038e+04 7217 0.9646 0.007005 0.03541 0.03519\n", " 13 0.54 2.087e+04 8436 0.9596 0.006922 0.04042 0.04011\n", " 14 0.58 2.174e+04 9842 0.9547 0.006782 0.04526 0.04513\n", " 15 0.62 2.446e+04 1.232e+04 0.9496 0.006394 0.05037 0.0502\n", " 16 0.66 2.577e+04 1.428e+04 0.9446 0.006229 0.05541 0.05527\n", " 17 0.7 2.662e+04 1.601e+04 0.9398 0.006129 0.06015 0.06029\n", " 18 0.74 2.361e+04 1.536e+04 0.9349 0.006508 0.06506 0.0652\n", " 19 0.78 2.285e+04 1.591e+04 0.9304 0.006615 0.06964 0.06995\n", " 20 0.82 1.933e+04 1.433e+04 0.9259 0.007192 0.07411 0.0745\n", " 21 0.86 2.055e+04 1.609e+04 0.9217 0.006976 0.07828 0.07879\n", " 22 0.9 2.038e+04 1.679e+04 0.9176 0.007004 0.08237 0.08279\n", " 23 0.94 2.36e+04 2.023e+04 0.9143 0.006509 0.08572 0.08645\n", " 24 0.98 2.149e+04 1.922e+04 0.9105 0.006822 0.08946 0.08973\n", " 25 1.02 2.046e+04 1.885e+04 0.9079 0.006991 0.09212 0.09261\n", " 26 1.06 2.049e+04 1.936e+04 0.9056 0.006985 0.09445 0.09504\n", " 27 1.1 1.631e+04 1.571e+04 0.9037 0.00783 0.09631 0.09701\n", " 28 1.14 1.72e+04 1.682e+04 0.9022 0.007624 0.09776 0.09849\n", " 29 1.18 1.586e+04 1.564e+04 0.9014 0.007941 0.09862 0.09948\n", " 30 1.22 1.716e+04 1.699e+04 0.901 0.007633 0.09902 0.09995\n", " 31 1.26 1.175e+04 1.162e+04 0.901 0.009227 0.09897 0.09991\n", " 32 1.3 1.196e+04 1.177e+04 0.9016 0.009143 0.0984 0.09936\n", " 33 1.34 9030 8783 0.9027 0.01052 0.09726 0.09829\n", " 34 1.38 1.052e+04 1.006e+04 0.9044 0.009748 0.09558 0.09673\n", " 35 1.42 8524 7967 0.9065 0.01083 0.09347 0.09469\n", " 36 1.46 8091 7360 0.909 0.01112 0.09096 0.09218\n", " 37 1.5 6648 5865 0.9118 0.01226 0.08823 0.08924\n", " 38 1.54 4921 4162 0.9154 0.01426 0.08458 0.0859\n", " 39 1.58 4600 3734 0.9188 0.01474 0.08118 0.08218\n", " 40 1.62 5665 4371 0.9228 0.01329 0.07716 0.07813\n", " 41 1.66 6462 4707 0.9272 0.01244 0.07284 0.0738\n", " 42 1.7 5763 3881 0.9326 0.01317 0.06735 0.06921\n", " 43 1.74 4660 2970 0.9363 0.01465 0.06374 0.06443\n", " 44 1.78 3583 2133 0.9405 0.01671 0.05952 0.0595\n", " 45 1.82 3369 1795 0.9467 0.01723 0.0533 0.05447\n", " 46 1.86 2767 1362 0.9508 0.01901 0.04922 0.0494\n", " 47 1.9 1207 529.5 0.9561 0.02879 0.04388 0.04433\n", " 48 1.94 2834 1108 0.9609 0.01879 0.03909 0.03932\n", " 49 1.98 2237 788 0.9648 0.02114 0.03522 0.03442\n", " 50 2.02 2056 629 0.9694 0.02205 0.03059 0.02968\n", " 51 2.06 1603 406.9 0.9746 0.02497 0.02538 0.02516\n", " 52 2.1 1493 325.4 0.9782 0.02588 0.02179 0.02088\n", " 53 2.14 1897 347.8 0.9817 0.02296 0.01833 0.01691\n", " 54 2.18 512.7 75.02 0.9854 0.04416 0.01463 0.01328\n", " 55 2.22 980.5 125 0.9873 0.03194 0.01274 0.01003\n", " 56 2.26 1006 100.7 0.99 0.03153 0.01001 0.007186\n", " 57 2.3 784.2 59.38 0.9924 0.03571 0.007572 0.004788\n", " 58 2.34 745.5 42.81 0.9943 0.03662 0.005742 0.002855\n", " 59 2.38 719.5 33.14 0.9954 0.03728 0.004606 0.001409\n", " 60 2.42 438.7 17.62 0.996 0.04775 0.004016 0.0004643\n", " 61 2.46 1076 40.31 0.9963 0.03049 0.003746 3.025e-05\n", " 62 2.5 421.5 18.44 0.9956 0.04871 0.004375 0.0001116\n", " 63 2.54 545.3 27.08 0.995 0.04282 0.004966 0.0007074\n", " 64 2.58 717.3 44.34 0.9938 0.03734 0.006181 0.001812\n", " 65 2.62 555.1 44.42 0.992 0.04244 0.008002 0.003413\n", " 66 2.66 441.6 43.59 0.9901 0.04759 0.009873 0.005495\n", " 67 2.7 202.4 26.27 0.987 0.07029 0.01298 0.008035\n", " 68 2.74 230.2 36.37 0.9842 0.06591 0.0158 0.01101\n", " 69 2.78 141.6 25.64 0.9819 0.08403 0.0181 0.01438\n", " 70 2.82 316.1 69.37 0.9781 0.05625 0.02195 0.01813\n", " 71 2.86 158.3 40.93 0.9741 0.07949 0.02586 0.0222\n", " 72 2.9 239.4 73.3 0.9694 0.06463 0.03062 0.02656\n", " 73 2.94 132.7 44.78 0.9663 0.0868 0.03374 0.03116\n", " 74 2.98 111.4 43.88 0.9606 0.09476 0.03941 0.03595\n", " 75 3.02 183.3 78.84 0.957 0.07386 0.043 0.04089\n", " 76 3.06 111.1 51.49 0.9536 0.09489 0.04636 0.04593\n", " 77 3.1 161.6 85 0.9474 0.07867 0.05261 0.051\n", " 78 3.14 91.11 51.8 0.9431 0.1048 0.05685 0.05607\n", " 79 3.18 163.3 97.15 0.9405 0.07825 0.05949 0.06107\n", " 80 3.22 64.97 41.94 0.9354 0.1241 0.06456 0.06596\n", " 81 3.26 122.3 83.82 0.9315 0.09041 0.06852 0.07068\n", " 82 3.3 93.82 68.36 0.9271 0.1032 0.07286 0.07519\n", " 83 3.34 136.9 104.3 0.9238 0.08546 0.07617 0.07944\n", " 84 3.38 195.4 153.5 0.9214 0.07155 0.07857 0.08339\n", " 85 3.42 109 87.96 0.9193 0.09579 0.08071 0.08699\n", " 86 3.46 65.4 55.15 0.9157 0.1237 0.08432 0.09021\n", " 87 3.5 37.29 32.33 0.9133 0.1637 0.0867 0.09302\n", " 88 3.54 80.11 70.54 0.9119 0.1117 0.08805 0.09538\n", " 89 3.58 113.9 101.3 0.911 0.0937 0.08898 0.09728\n", " 90 3.62 59.96 54.06 0.9098 0.1291 0.09016 0.09868\n", " 91 3.66 19.93 18.03 0.9095 0.224 0.09046 0.09959\n", " 92 3.7 20.89 18.9 0.9095 0.2188 0.0905 0.09998\n", " 93 3.74 4.295 3.677 0.9144 0.4825 0.08562 0.09986\n", " 94 3.78 46.93 41.18 0.9123 0.146 0.08775 0.09922\n", " 95 3.82 86.48 75.93 0.9122 0.1075 0.0878 0.09808\n", " 96 3.86 46.69 40.19 0.9139 0.1463 0.08607 0.09644\n", " 97 3.9 26.84 22.4 0.9165 0.193 0.08345 0.09432\n", " 98 3.94 62.32 50.84 0.9184 0.1267 0.08159 0.09175\n", " 99 3.98 30.15 23.75 0.9212 0.1821 0.07877 0.08874\n", "100 4.02 31.25 24.12 0.9228 0.1789 0.07716 0.08533\n", "101 4.06 19.22 14.19 0.9262 0.2281 0.07383 0.08156\n", "102 4.1 9.354 6.371 0.9319 0.327 0.06811 0.07747\n", "103 4.14 1.905 1.239 0.935 0.7246 0.06504 0.07309\n", "104 4.18 35.17 22.23 0.9368 0.1686 0.06322 0.06847\n", "105 4.22 36 21.16 0.9412 0.1667 0.05879 0.06366\n", "106 4.26 11.19 6.198 0.9446 0.2989 0.05538 0.05871\n", "107 4.3 8.742 4.474 0.9488 0.3382 0.05118 0.05368\n", "108 4.34 7.641 3.743 0.951 0.3618 0.04898 0.0486\n", "109 4.38 25.41 11.14 0.9561 0.1984 0.04386 0.04354\n", "110 4.42 27.43 11.35 0.9586 0.1909 0.0414 0.03854\n", "111 4.46 9.079 3.426 0.9623 0.3319 0.03774 0.03366\n", "112 4.5 0.2087 0.07581 0.2011 1 0.007581 0.02896\n", "113 4.54 22.14 7.139 0.9678 0.2125 0.03224 0.02446\n", "114 4.58 2.152 0.6842 0.9682 0.6817 0.03179 0.02024\n", "115 4.62 5.036 1.458 0.9711 0.4456 0.02895 0.01631\n", "116 4.66 31.57 7.739 0.9755 0.178 0.02451 0.01274\n", "117 4.7 7.329 1.802 0.9754 0.3694 0.02459 0.009551\n", "118 4.74 3.533 0.7965 0.9775 0.532 0.02254 0.006778\n", "119 4.78 0.2103 0.05084 0.2052 1 0.005084 0.004452\n", "120 4.82 18 3.398 0.9811 0.2357 0.01888 0.002595\n", "121 4.86 6.921 1.232 0.9822 0.3801 0.0178 0.001227\n", "122 4.9 3.444 0.6831 0.9802 0.5389 0.01984 0.0003618\n", "123 4.94 14.82 2.616 0.9824 0.2597 0.01765 8.808e-06\n", "124 4.98 10.14 2.138 0.9789 0.3141 0.02109 0.0001714\n", "125 5.02 0.9335 0.2223 0.9113 1 0.02223 0.0008479\n", "126 5.06 6.348 1.356 0.9786 0.3969 0.02136 0.002031\n", "127 5.1 3.209 0.7468 0.9767 0.5582 0.02327 0.00371\n", "128 5.14 3.229 0.8033 0.9751 0.5565 0.02487 0.005865\n", "129 5.18 4.162 1.132 0.9728 0.4902 0.02721 0.008476\n", "130 5.22 0 0 0 1 0 0.01152\n", "131 5.26 2.849 0.9148 0.9679 0.5925 0.03211 0.01495\n", "132 5.3 6.447 2.111 0.9673 0.3938 0.03275 0.01875\n", "133 5.34 3.217 1.113 0.9654 0.5575 0.03461 0.02287\n", "134 5.38 0.9703 0.3801 0.9323 1 0.03801 0.02727\n", "135 5.42 3.221 1.296 0.9598 0.5572 0.04023 0.0319\n", "136 5.46 0 0 0 1 0 0.03672\n", "137 5.5 0 0 0 1 0 0.04168\n", "138 5.54 3.229 1.531 0.9526 0.5565 0.0474 0.04672\n", "139 5.58 3.126 1.616 0.9483 0.5656 0.05171 0.0518\n", "140 5.62 3.225 1.779 0.9448 0.5568 0.05516 0.05686\n", "141 5.66 0.9728 0.535 0.9193 1 0.0535 0.06185\n", "142 5.7 0 0 0 1 0 0.06671\n", "143 5.74 0 0 0 1 0 0.07141\n", "144 5.78 0.9691 0.5921 0.9099 1 0.05921 0.07588\n", "145 5.82 0.9716 0.6185 0.9097 1 0.06185 0.08009\n", "146 5.86 0.9655 0.6205 0.9035 1 0.06205 0.08398\n", "147 5.9 1.874 1.243 0.9337 0.7305 0.06634 0.08753\n", "148 5.94 0 0 0 1 0 0.09068\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", "nyrs = 5.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 = 150\n", "# nle = 30\n", "dle_max = 6.\n", "\n", "# print out run parameters\n", "with open(\"current_istope_event_calc_run_CHANDLERW3.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_CHANDLERW3.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", "print 'nev: ',nev\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", " #print 'iev: ',iev\n", " line = contents[iev+1].strip() # changed from -1 to +1\n", " #print line\n", " # enu,levt,n_nuebar_IBD,n_nuebar\n", " fields = line.split()\n", " #print fields\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", " #print dat_array[0,iev],dat_array[1,iev],dat_array[3,iev]\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(10,71):\n", "for idm2 in range(40,41):\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 True:\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", " dle = dle_max/nle\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 }