{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the python equivalent of isotope_event_calculator.f written by Mike Shaevitz. It was adapted for python by William Seligman. \n",
"\n",
"_Markdown cells in italics indicate a note added by Seligman during his conversion process. First note: The odds that I've made an arithmetic error in the code conversion is probably 100%._\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": "markdown",
"metadata": {},
"source": [
"_First, define a function call for use by the main program. I took the liberty of splitting up the polynomial into separate lines to make sure I copied the terms over correctly._"
]
},
{
"cell_type": "code",
"execution_count": 287,
"metadata": {},
"outputs": [],
"source": [
"# IBD cross section in cm^2\n",
"def IBD_xsec(eee):\n",
" xsec = (\n",
" 1.08497e-06 \n",
" - 0.00971813*eee\n",
" + 0.00484432*eee**2\n",
" + 0.000521442*eee**3 \n",
" - 2.81264e-05*eee**4\n",
" + 5.62992e-07*eee**5 \n",
" - 3.90173e-09*eee**6\n",
" )\n",
" return max(0.0,xsec)*1.e-41"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Initialize constants. Note that floating-point variables in Python are double-precision (\"float64\") by default._"
]
},
{
"cell_type": "code",
"execution_count": 288,
"metadata": {},
"outputs": [],
"source": [
"nenu = 100\n",
"nl = 100\n",
"ny = 100\n",
"\n",
"nebins = 50\n",
"nlbins = 50\n",
"emax_ib = 20.\n",
"\n",
"shimi_emax = 15. # only use values from file less than this\n",
"\n",
"# inputs to calculation\n",
"sinsqthw = 0.238"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Python can only write strings to output files, so you have to convert any numbers before you write them. Also, when writing to a file, you have explicitly include \"\\n\" (newline) at the end of a line unless you want the next .write() to add to the end of the previous output line._"
]
},
{
"cell_type": "code",
"execution_count": 289,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"10000 1000000 0.238\n"
]
}
],
"source": [
"eventfile = open(\"nuebar_event_array.dat\",\"w\")\n",
"eventfile.write( str(nenu*nl) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"evisfile = open(\"nuebar_evis_array.dat\",\"w\")\n",
"evisfile.write( str(nenu*nl*ny) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"calcrunfile = open(\"current_istope_event_calc_run.txt\",\"w\")\n",
"line = str(nenu*nl) + \" \" + str(nenu*nl*ny) + \" \" + str(sinsqthw)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line"
]
},
{
"cell_type": "code",
"execution_count": 290,
"metadata": {},
"outputs": [],
"source": [
"# 5mA of alphas is:\n",
"# (5e-3 C/s)*(1./2*1.6e-19 n_alpha/C)=1.56e16/s\n",
"# (1.56e16/s)*365*24*3600=4.93e23 alphas/yr\n",
"# 5ma of H2+ (or 10ms of protons):\n",
"# (5e-3 C/s)*(2./1.6e-19 H2+/C)=6.25e16p/s\n",
"# protons/yr = 6.25e16p/s*365*24*3600 = \n",
"protons_per_year = 1.97e24 # 10ma_protons @ 100%df\n",
"\n",
"# isotope production by alphas: include 8Li and 9Li 240 MeV alphas from Adriana\n",
"# Rate(8Li+9Li)/alpha) = (71827+5643)/1.e7 = 0.00775\n",
"# isotopes_per_proton = 8.9e-4 # 12N see Estimates of Interaction and Production Rates.xls\n",
"# isotope production by protons 8Li is 0.0014/proton\n",
"# isotopes_per_proton = 0.0014 # Adriana email 2/24/2012 20cm thick target proton production\n",
"# isotopes_per_proton = 0.00152 # Adriana update extended sleeve geometry 2/29/12 60 MeV\n",
"# isotopes_per_proton = 0.0045 # x3 Adriana update extended sleeve geometry 2/29/12 60 MeV\n",
"# isotopes_per_proton = 0.0023 # x1.5 for pure 11B sleeve Adriana 2/29/12 60 MeV\n",
"# isotopes_per_proton = (13495.+3827.+281221.)/1.e7 # 0.03 Adriana email 3/12/12 60 MeV\n",
"# isotopes_per_proton = (13677. + 121.)/1.e7 # Back to 8Li only 60 MeV\n",
"# isotopes_per_proton = 4.*(13677. + 121.)/1.e7 # x4 8Li only rate\n",
"isotopes_per_proton = 145600/1.e7 # Adriana email 4/30/12 Pure 9Be + 99.99 7Li\n",
"\n",
"isotope_br = 1.0 \n",
"\n",
"# nyrs =3.82e21/(protons_per_year*isotopes_per_proton*isotope_br)\n",
"# nyrs = 3.0*0.8 # 3yrs at 80% DF\n",
"# nyrs = 1.0*0.8 # 1yrs at 80% DF\n",
"# 80% DF includes 90% DF plus 92% detection efficiency (in the PRL)\n",
"nyrs = 5.0*0.8 # 5yrs at 80% DF\n",
"\n",
"total_flux=nyrs*protons_per_year*isotopes_per_proton*isotope_br"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Formatting output in Python is not the same as in FORTRAN. The floating-point formatting specifications contains the same characters but organized differently; FORTRAN's _ f8.3 _becomes_ {:8.3f}. For integers, i5 becomes {:5d}. \n",
"\n",
"_Instead of FORTRAN's FORMAT statement, in Python we create a string with the format then use a format() function on the string to supply the variables to be formatted. Example:_\n",
"\n",
" write (6,20) myvalue\n",
" 20 format(\"My value is \",f8.3\")\n",
" \n",
"_becomes in Python_\n",
"\n",
" print (\"My value is {:8.3f}\".format(myvalue))\n",
"\n",
"_Because we're writing the same thing to two different files, I'm using Python's string formatting to create a text string (\"line\") then writing that line to both files._"
]
},
{
"cell_type": "code",
"execution_count": 291,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total flux: 1.147e+23 nyrs(with DF): 4.00 protons/yr: 1.97e+24\n",
"isotopes/proton: 0.01456 isotope_br: 1.00\n",
"beta-decay endpoint = 15.000 MeV\n",
"mass: 897.0 radius: 6.50 shield: 3.2 distance: 12.50\n"
]
}
],
"source": [
"from math import pi\n",
"\n",
"line = \"total flux: {:12.4g} nyrs(with DF): {:5.2f} protons/yr: {:12.4g}\".format(total_flux,nyrs,protons_per_year)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line\n",
"\n",
"line = \"isotopes/proton: {:9.5f} isotope_br: {:5.2f}\".format(isotopes_per_proton,isotope_br)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line\n",
"\n",
"# e_endpoint = 16.3 # N12 beta-decay endpoint (MeV)\n",
"# e_endpoint = 16.1 # 8Li beta-decay endpoint (MeV) (used before 3/1/12)\n",
"# e_endpoint = 12.6 # 8Li beta-decay endpoint (MeV) from fit to Shimizu electron spectra\n",
"e_endpoint = 15.0 # emax for Shimizu nuebar flux file\n",
"emax_ib = 15.0\n",
"\n",
"line = \"beta-decay endpoint = {:8.3f} MeV\".format(e_endpoint)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line\n",
"\n",
"# detector_mass = 1000. # metric tons (used before 3/1/12)\n",
"# detector_mass = 800. # metric tons corresponds to a 6.1m detector radius\n",
"# detector_mass = 1000. # metric tons - use full volume of KamLand\n",
"detector_mass = 897. # metric tons - use full volume of KamLand\n",
"\n",
"# detector_density = 0.85 # g/cm3\n",
"detector_density = 0.78 # g/cm3\n",
"\n",
"detector_volume=(detector_mass*1000.)/(detector_density*1000.)\n",
"detector_radius = (.75*detector_mass*1000./(pi*detector_density*1000.))**(1./3.)\n",
"# lshield = 10. # lenght of shielding (m)\n",
"# lshield = 5. # lenght of shielding (m)\n",
"# d_source_detector = lshield+detector_radius # distance from source to detector (m)\n",
"# d_source_detector = 16.5 # distance from source to detector (m) <== Conservative\n",
"# d_source_detector = 16.0 # distance from source to detector (m) <== Agressive <= nominal\n",
"d_source_detector = 12.5 # distance from source to detector (m) <== New Option 2 possibility\n",
"# lshield = d_source_detector-13.0 # 9m containment vessel plus 1m to rock.\n",
"lshield = d_source_detector-9.25 # Det_Rad+buffer+PMT_sphere+concrete = 6.5+2.2+0.55=9.25\n",
"\n",
"line = \"mass: {:8.1f} radius: {:6.2f} shield: {:6.1f} distance: {:6.2f}\".format(detector_mass,detector_radius,lshield,d_source_detector)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_These are Python \"lists\" as opposed to numpy arrays. It's easier to append values to the end of lists. Start with empty lists._\n"
]
},
{
"cell_type": "code",
"execution_count": 292,
"metadata": {},
"outputs": [],
"source": [
"tmp_e = []\n",
"tmp_flux = []"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Read the shmizu file. All Python I/O is based on strings, so we have to convert the lines in the file into numeric values. Typically, once we are done with a file, we must close it. (We haven't close the files we've opened above because we've got more to write.)\n",
"\n",
"_Here we're using the \"with\" method of reading the file, because after reading an entire file using \"with\" the file is automatically closed._\n",
"\n",
"_Reminder: loops, ifs, and other control-flow statements are indicated by indents in Python._"
]
},
{
"cell_type": "code",
"execution_count": 293,
"metadata": {},
"outputs": [],
"source": [
"# read in Shimizu file\n",
"with open(\"antinus_8li_shimizu.dat\",\"r\") as f:\n",
" contents = f.readlines() # Read the file into a list of lines.\n",
" \n",
"nshimi_count = len(contents) # the total number of lines in the shimizu file\n",
"\n",
"for line in contents:\n",
" line = line.strip() # Remove end-of-line marker\n",
" e, flux = line.split() # Extract the two non-blank fields from the string\n",
" energy = float(e) # Convert the string into a number\n",
" if ( energy < shimi_emax ):\n",
" tmp_e.append( energy ) # Append the values to the end of the list\n",
" tmp_flux.append( float(flux)/100. )\n",
" \n",
"nshimi_file = len(tmp_e) # The number of entries in the shimizu file that pass the cut"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Note that I've used arrays from the **numpy** package. There are other ways to handle 1-dimensional arrays in Python, but numpy arrays can be used as arguments to ROOT functions, which may be handy in the future._"
]
},
{
"cell_type": "code",
"execution_count": 294,
"metadata": {},
"outputs": [],
"source": [
"# Use numpy arrays\n",
"import numpy as np\n",
"\n",
"# Create an 1-dimensional array of \"nenu\" elements all set to zero.\n",
"shimi_e = np.zeros((nenu))\n",
"shimi_flux = np.zeros((nenu))\n",
"\n",
"nevts_nue_e = np.zeros((nebins))\n",
"nevts_nuebar_e = np.zeros((nebins))\n",
"nevts_nuebar_IBD = np.zeros((nebins))\n",
"nvis_nue_e = np.zeros((nebins))\n",
"nvis_nuebar_e = np.zeros((nebins))\n",
"nevts_nuebar_IBD_levt = np.zeros((nlbins))"
]
},
{
"cell_type": "code",
"execution_count": 295,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Shimizu flux dist\n",
" 0 0.075 2.41757e-05\n",
" 1 0.225 0.00013198\n",
" 2 0.375 0.000338553\n",
" 3 0.525 0.000634859\n",
" 4 0.675 0.00100964\n",
" 5 0.825 0.00145464\n",
" 6 0.975 0.00196037\n",
" 7 1.125 0.00251963\n",
" 8 1.275 0.00312429\n",
" 9 1.425 0.00376794\n",
" 10 1.575 0.00444356\n",
" 11 1.725 0.0051454\n",
" 12 1.875 0.00586738\n",
" 13 2.025 0.00660435\n",
" 14 2.175 0.00735103\n",
" 15 2.325 0.00810278\n",
" 16 2.475 0.00885509\n",
" 17 2.625 0.00960376\n",
" 18 2.775 0.0103449\n",
" 19 2.925 0.0110749\n",
" 20 3.075 0.0117903\n",
" 21 3.225 0.0124879\n",
" 22 3.375 0.0131649\n",
" 23 3.525 0.0138184\n",
" 24 3.675 0.0144462\n",
" 25 3.825 0.0150457\n",
" 26 3.975 0.0156151\n",
" 27 4.125 0.0161522\n",
" 28 4.275 0.0166557\n",
" 29 4.425 0.0171237\n",
" 30 4.575 0.0175552\n",
" 31 4.725 0.0179487\n",
" 32 4.875 0.0183037\n",
" 33 5.025 0.0186188\n",
" 34 5.175 0.0188939\n",
" 35 5.325 0.019128\n",
" 36 5.475 0.0193212\n",
" 37 5.625 0.0194728\n",
" 38 5.775 0.0195832\n",
" 39 5.925 0.0196521\n",
" 40 6.075 0.0196801\n",
" 41 6.225 0.019667\n",
" 42 6.375 0.0196138\n",
" 43 6.525 0.0195207\n",
" 44 6.675 0.0193888\n",
" 45 6.825 0.0192185\n",
" 46 6.975 0.0190111\n",
" 47 7.125 0.0187673\n",
" 48 7.275 0.0184886\n",
" 49 7.425 0.018176\n",
" 50 7.575 0.0178311\n",
" 51 7.725 0.017455\n",
" 52 7.875 0.0170496\n",
" 53 8.025 0.0166163\n",
" 54 8.175 0.016157\n",
" 55 8.325 0.0156732\n",
" 56 8.475 0.0151671\n",
" 57 8.625 0.0146404\n",
" 58 8.775 0.0140953\n",
" 59 8.925 0.0135336\n",
" 60 9.075 0.0129577\n",
" 61 9.225 0.0123697\n",
" 62 9.375 0.0117718\n",
" 63 9.525 0.0111663\n",
" 64 9.675 0.0105556\n",
" 65 9.825 0.00994185\n",
" 66 9.975 0.0093276\n",
" 67 10.125 0.00871522\n",
" 68 10.275 0.00810708\n",
" 69 10.425 0.00750566\n",
" 70 10.575 0.00691331\n",
" 71 10.725 0.00633253\n",
" 72 10.875 0.00576557\n",
" 73 11.025 0.00521496\n",
" 74 11.175 0.0046828\n",
" 75 11.325 0.00417157\n",
" 76 11.475 0.0036832\n",
" 77 11.625 0.00322004\n",
" 78 11.775 0.00278377\n",
" 79 11.925 0.00237656\n",
" 80 12.075 0.00199972\n",
" 81 12.225 0.00165516\n",
" 82 12.375 0.00134376\n",
" 83 12.525 0.00106701\n",
" 84 12.675 0.000825269\n",
" 85 12.825 0.000619504\n",
" 86 12.975 0.00044934\n",
" 87 13.125 0.000314506\n",
" 88 13.275 0.000212366\n",
" 89 13.425 0.000139161\n",
" 90 13.575 8.88055e-05\n",
" 91 13.725 5.55627e-05\n",
" 92 13.875 3.40266e-05\n",
" 93 14.025 2.0433e-05\n",
" 94 14.175 1.19419e-05\n",
" 95 14.325 6.78146e-06\n",
" 96 14.475 3.69641e-06\n",
" 97 14.625 1.92459e-06\n",
" 98 14.775 9.38851e-07\n"
]
}
],
"source": [
"# rebin into nenu energy bins\n",
"nshimi = nshimi_file/nenu\n",
"print 'Shimizu flux dist'\n",
"for ienu in range (0,nenu-1):\n",
" for ib in range (ienu*nshimi, (ienu+1)*nshimi-1):\n",
" shimi_e[ienu] = shimi_e[ienu] + tmp_e[ib]*tmp_flux[ib]\n",
" shimi_flux[ienu] = shimi_flux[ienu] + tmp_flux[ib]\n",
" if shimi_flux[ienu] > 0. :\n",
" shimi_e[ienu] = shimi_e[ienu]/shimi_flux[ienu]\n",
" else:\n",
" ib = int((ienu-0.5)*nshimi) - 1\n",
" shimi_e[ienu] = tmp_e[ib]+0.005\n",
" shimi_e[ienu] = (shimi_emax/nenu)*(ienu+0.5)\n",
" print \"{:3d}{:14.6g}{:14.6}\".format(ienu,shimi_e[ienu],shimi_flux[ienu])"
]
},
{
"cell_type": "code",
"execution_count": 296,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"enu,wt_spectrum,flux0,IBD_xsec(enu)\n",
" 0.075 2.418e-05 2.774e+18 0\n",
" 0.225 0.000132 1.514e+19 0\n",
" 0.375 0.0003386 3.884e+19 0\n",
" 0.525 0.0006349 7.284e+19 0\n",
" 0.675 0.00101 1.158e+20 0\n",
" 0.825 0.001455 1.669e+20 0\n",
" 0.975 0.00196 2.249e+20 0\n",
" 1.125 0.00252 2.891e+20 0\n",
" 1.275 0.003124 3.585e+20 0\n",
" 1.425 0.003768 4.323e+20 0\n",
" 1.575 0.004444 5.098e+20 0\n",
" 1.725 0.005145 5.903e+20 8.818e-46\n",
" 1.875 0.005867 6.732e+20 1.913e-44\n",
" 2.025 0.006604 7.577e+20 4.062e-44\n",
" 2.175 0.007351 8.434e+20 6.544e-44\n",
" 2.325 0.008103 9.297e+20 9.362e-44\n",
" 2.475 0.008855 1.016e+21 1.252e-43\n",
" 2.625 0.009604 1.102e+21 1.604e-43\n",
" 2.775 0.01034 1.187e+21 1.99e-43\n",
" 2.925 0.01107 1.271e+21 2.413e-43\n",
" 3.075 0.01179 1.353e+21 2.872e-43\n",
" 3.225 0.01249 1.433e+21 3.368e-43\n",
" 3.375 0.01316 1.51e+21 3.902e-43\n",
" 3.525 0.01382 1.585e+21 4.473e-43\n",
" 3.675 0.01445 1.657e+21 5.083e-43\n",
" 3.825 0.01505 1.726e+21 5.731e-43\n",
" 3.975 0.01562 1.792e+21 6.419e-43\n",
" 4.125 0.01615 1.853e+21 7.145e-43\n",
" 4.275 0.01666 1.911e+21 7.911e-43\n",
" 4.425 0.01712 1.965e+21 8.718e-43\n",
" 4.575 0.01756 2.014e+21 9.564e-43\n",
" 4.725 0.01795 2.059e+21 1.045e-42\n",
" 4.875 0.0183 2.1e+21 1.138e-42\n",
" 5.025 0.01862 2.136e+21 1.235e-42\n",
" 5.175 0.01889 2.168e+21 1.336e-42\n",
" 5.325 0.01913 2.195e+21 1.441e-42\n",
" 5.475 0.01932 2.217e+21 1.55e-42\n",
" 5.625 0.01947 2.234e+21 1.663e-42\n",
" 5.775 0.01958 2.247e+21 1.781e-42\n",
" 5.925 0.01965 2.255e+21 1.902e-42\n",
" 6.075 0.01968 2.258e+21 2.028e-42\n",
" 6.225 0.01967 2.256e+21 2.158e-42\n",
" 6.375 0.01961 2.25e+21 2.292e-42\n",
" 6.525 0.01952 2.24e+21 2.431e-42\n",
" 6.675 0.01939 2.225e+21 2.573e-42\n",
" 6.825 0.01922 2.205e+21 2.72e-42\n",
" 6.975 0.01901 2.181e+21 2.871e-42\n",
" 7.125 0.01877 2.153e+21 3.026e-42\n",
" 7.275 0.01849 2.121e+21 3.186e-42\n",
" 7.425 0.01818 2.085e+21 3.349e-42\n",
" 7.575 0.01783 2.046e+21 3.517e-42\n",
" 7.725 0.01745 2.003e+21 3.689e-42\n",
" 7.875 0.01705 1.956e+21 3.865e-42\n",
" 8.025 0.01662 1.906e+21 4.045e-42\n",
" 8.175 0.01616 1.854e+21 4.23e-42\n",
" 8.325 0.01567 1.798e+21 4.418e-42\n",
" 8.475 0.01517 1.74e+21 4.611e-42\n",
" 8.625 0.01464 1.68e+21 4.807e-42\n",
" 8.775 0.0141 1.617e+21 5.008e-42\n",
" 8.925 0.01353 1.553e+21 5.213e-42\n",
" 9.075 0.01296 1.487e+21 5.422e-42\n",
" 9.225 0.01237 1.419e+21 5.635e-42\n",
" 9.375 0.01177 1.351e+21 5.852e-42\n",
" 9.525 0.01117 1.281e+21 6.073e-42\n",
" 9.675 0.01056 1.211e+21 6.298e-42\n",
" 9.825 0.009942 1.141e+21 6.526e-42\n",
" 9.975 0.009328 1.07e+21 6.759e-42\n",
" 10.12 0.008715 9.999e+20 6.996e-42\n",
" 10.28 0.008107 9.301e+20 7.236e-42\n",
" 10.42 0.007506 8.611e+20 7.481e-42\n",
" 10.57 0.006913 7.932e+20 7.729e-42\n",
" 10.72 0.006333 7.265e+20 7.981e-42\n",
" 10.88 0.005766 6.615e+20 8.237e-42\n",
" 11.03 0.005215 5.983e+20 8.496e-42\n",
" 11.17 0.004683 5.373e+20 8.759e-42\n",
" 11.32 0.004172 4.786e+20 9.026e-42\n",
" 11.47 0.003683 4.226e+20 9.297e-42\n",
" 11.62 0.00322 3.694e+20 9.571e-42\n",
" 11.78 0.002784 3.194e+20 9.849e-42\n",
" 11.92 0.002377 2.727e+20 1.013e-41\n",
" 12.07 0.002 2.294e+20 1.042e-41\n",
" 12.22 0.001655 1.899e+20 1.07e-41\n",
" 12.38 0.001344 1.542e+20 1.1e-41\n",
" 12.53 0.001067 1.224e+20 1.129e-41\n",
" 12.67 0.0008253 9.469e+19 1.159e-41\n",
" 12.82 0.0006195 7.108e+19 1.189e-41\n",
" 12.97 0.0004493 5.155e+19 1.22e-41\n",
" 13.12 0.0003145 3.608e+19 1.251e-41\n",
" 13.28 0.0002124 2.437e+19 1.282e-41\n",
" 13.42 0.0001392 1.597e+19 1.313e-41\n",
" 13.57 8.881e-05 1.019e+19 1.345e-41\n",
" 13.72 5.556e-05 6.375e+18 1.377e-41\n",
" 13.88 3.403e-05 3.904e+18 1.41e-41\n",
" 14.03 2.043e-05 2.344e+18 1.443e-41\n",
" 14.17 1.194e-05 1.37e+18 1.476e-41\n",
" 14.32 6.781e-06 7.781e+17 1.509e-41\n",
" 14.47 3.696e-06 4.241e+17 1.543e-41\n",
" 14.62 1.925e-06 2.208e+17 1.577e-41\n",
" 14.77 9.389e-07 1.077e+17 1.611e-41\n"
]
}
],
"source": [
"from math import sqrt # In python, this is not built-in.\n",
"\n",
"# Loop over beta-decay spectrum\n",
"# denu = e_endpoint/nenu\n",
"print 'enu,wt_spectrum,flux0,IBD_xsec(enu)'\n",
"for ienu in range (0,nenu-1):\n",
" enu = shimi_e[ienu]\n",
" wt_spectrum = shimi_flux[ienu]\n",
" # enu = (ienu-0.5)*denu\n",
" # wt_spectrum = max(0.e0,enu**2*(e_endpoint-enu)**2*denu)\n",
" # wt_spectrum = wt_spectrum/((e_endpoint**5)/30.) # normalize to 1.0\n",
" flux0 = protons_per_year*isotopes_per_proton*nyrs*isotope_br*wt_spectrum \n",
" # print 'enu,wt_spectrum,flux0: ',enu,wt_spectrum,flux0\n",
"\n",
" # IBD xsec for antineutrinos\n",
" ee = max(0.511,enu - 1.29332)\n",
" pe = sqrt(ee**2-0.511**2)\n",
" # xsec (10**-43 cm2 from Vogel/Beacom scale to xsec1 from xsec0)\n",
" # need to put in factor (6_protons/8_electrons) to be consistent with nu-electron xsec\n",
" # ### Mistake: should be factor (2_free_protons/8_electrons) to be consistent with nu-electron xsec\n",
" # xsecIBD = 10.*0.0952*ee*pe*1.e-43*(6./8.)\n",
" # xsecIBD = 10.*0.0952*ee*pe*1.e-43*2.0 # <= correct 2 free protons per CH2\n",
" xsecIBD = IBD_xsec(enu)*2.0 # use parameterization with correct 2 free protons per CH2\n",
" if enu <= 1.29332: \n",
" xsecIBD = 0.0\n",
"\n",
" print \"{:12.4g}{:12.4g}{:12.4g}{:12.4g}\".format(enu,wt_spectrum,flux0,IBD_xsec(enu))\n",
"\n",
" # nu-e scattering xsecs\n",
" eminus = -0.5-sinsqthw\n",
" eplus = -sinsqthw\n",
" # sig0 = 2Gmu^2*me*enu/(pi) in cm^2\n",
" sig0 = 1.5e-44*enu # enu in MeV\n",
" xsec_nue_e = sig0*(eminus**2+eplus**2/3.0)*8.0 # 8 electrons per CH2\n",
" xsec_nuebar_e = sig0*(eplus**2+eminus**2/3.0)*8.0 # 8 electrons per CH2\n",
" # print 'enu,xsecs - nue-e,nuebar-e,IBD: ',enu,xsec_nue_e,xsec_nuebar_e,xsecIBD \n",
"\n",
" # integrate over detector volume\n",
" dl = 2.*detector_radius/nl\n",
" for il in range (1,nl):\n",
" levt = dl*(il-0.5)+d_source_detector-detector_radius\n",
" xevt = (levt**2-detector_radius**2+d_source_detector**2)/(2.*d_source_detector)\n",
" yevt = sqrt(levt**2-xevt**2)\n",
" costh = sqrt(1.-(yevt/levt)**2)\n",
"\n",
" # calculate the flux into this solid angle\n",
" flux = flux0*0.5*(1-costh)\n",
" # print dl,levt,xevt,yevt,0.5*(1.-costh),flux\n",
"\n",
" # calculate the events/xsec for this dl \n",
" # (CH2: 8 electrons/CH2, 14g/mole) Use grams and cm units.\n",
" n_per_xsec = detector_density*(100.*dl)*(1./14.)*6.02e23*flux\n",
" # print 'n_per_xsec: ',n_per_xsec\n",
" n_nue = n_per_xsec*xsec_nue_e\n",
" n_nuebar = n_per_xsec*xsec_nuebar_e\n",
" n_nuebar_IBD = n_per_xsec*xsecIBD\n",
" # print 'n_nue,n_nuebar: ',n_nue,n_nuebar,n_nuebar_IBD\n",
"\n",
" # bin events according to neutrino energy\n",
" de_ib = emax_ib/nebins\n",
" ib = int(enu/de_ib)\n",
" nevts_nue_e[ib] = nevts_nue_e[ib] + n_nue\n",
" nevts_nuebar_e[ib] = nevts_nuebar_e[ib] + n_nuebar\n",
" nevts_nuebar_IBD[ib] = nevts_nuebar_IBD[ib] + n_nuebar_IBD\n",
"\n",
" # bin events according to distance from source\n",
" dl_ib = (2.*detector_radius)/nlbins\n",
" ib = int((levt-(d_source_detector-detector_radius))/dl_ib)\n",
" nevts_nuebar_IBD_levt[ib] = nevts_nuebar_IBD_levt[ib] + n_nuebar_IBD\n",
"\n",
" # write out event array with etrue\n",
" line = \"{:12.4g}{:12.4g}\".format(enu,levt)\n",
" for i in range (0,n_nuebar_IBD.size-1):\n",
" line += \"{:12.4g}\".format(n_nuebar_IBD[i])\n",
" for i in range (0,n_nuebar.size-1):\n",
" line += \"{:12.4g}\".format(n_nuebar[i])\n",
" eventfile.write(line + \"\\n\")\n",
"\n",
" # electron (visible) energy found from the y distribution\n",
" dy = 1.0/ny\n",
" for iy in range (1,ny):\n",
" y = (iy-0.5)*dy\n",
" e_elect = y*abs(enu)*1.00\n",
" norm = eminus**2+eplus**2/3.0\n",
" wtnue = (eminus**2+eplus**2*(1-y)**2)*dy*n_nue/norm\n",
" norm = eplus**2+eminus**2/3.0\n",
" wtnuebar = (eplus**2+eminus**2*(1-y)**2)*dy*n_nuebar/norm\n",
" # bin events according the outgoing electron energy\n",
" de_ib = emax_ib/nebins\n",
" ib = int(1. + e_elect/de_ib)\n",
" nvis_nue_e[ib] = nvis_nue_e[ib]+ wtnue\n",
" nvis_nuebar_e[ib] = nvis_nuebar_e[ib]+ wtnuebar\n",
" # write out event array with evis\n",
" line = \"{:12.4g}{:12.4g}{:12.4g}{:12.4g}{:12.4g}\\n\".format(enu,levt,e_elect,wtnue,wtnuebar)\n",
" evisfile.write(line + \"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Do a bit of column formatting for the final table, to show one way to get text into a fixed-width field. Note the continuation lines (a trailing \"\\\" with the following line indented)._"
]
},
{
"cell_type": "code",
"execution_count": 297,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" enu(MeV) n_nue_e n_nuebar_e n_nue_e_evisn_nuebar_e_evis n_nuebar_IBD\n",
" 0.15 0.51 0.22 0.00 0.00 0.00\n",
" 0.45 7.45 3.15 2.67 2.54 0.00\n",
" 0.75 30.44 12.86 2.66 2.29 0.00\n",
" 1.05 76.77 32.45 2.58 2.01 0.00\n",
" 1.35 151.29 63.95 2.59 1.83 0.00\n",
" 1.65 256.78 108.54 2.57 1.64 2.17\n",
" 1.95 394.29 166.66 2.50 1.45 182.06\n",
" 2.25 563.37 238.13 2.42 1.27 593.06\n",
" 2.55 762.31 322.22 2.42 1.16 1267.42\n",
" 2.85 988.37 417.77 2.34 1.02 2263.56\n",
" 3.15 1237.92 523.25 2.30 0.92 3632.52\n",
" 3.45 1506.65 636.83 2.22 0.81 5414.91\n",
" 3.75 1789.70 756.47 2.15 0.72 7638.56\n",
" 4.05 2081.81 879.94 2.04 0.63 10316.51\n",
" 4.35 2377.47 1004.91 1.95 0.56 13445.71\n",
" 4.65 2671.02 1128.99 1.87 0.49 17006.15\n",
" 4.95 2956.81 1249.79 1.79 0.44 20960.60\n",
" 5.25 3229.25 1364.95 1.66 0.38 25254.89\n",
" 5.55 3482.98 1472.19 1.63 0.35 29818.64\n",
" 5.85 3712.90 1569.37 1.41 0.28 34566.61\n",
" 6.15 3914.31 1654.51 1.35 0.26 39400.44\n",
" 6.45 4082.99 1725.81 1.27 0.23 44210.83\n",
" 6.75 4215.23 1781.70 1.16 0.20 48880.25\n",
" 7.05 4307.97 1820.90 1.02 0.17 53285.83\n",
" 7.35 4358.80 1842.39 0.94 0.15 57302.78\n",
" 7.65 4366.06 1845.45 0.86 0.13 60807.95\n",
" 7.95 4328.87 1829.74 0.77 0.11 63683.63\n",
" 8.25 4247.20 1795.22 0.66 0.09 65821.70\n",
" 8.55 4121.87 1742.24 0.57 0.08 67127.52\n",
" 8.85 3954.59 1671.53 0.50 0.06 67524.31\n",
" 9.15 3748.00 1584.21 0.43 0.05 66957.09\n",
" 9.45 3505.65 1481.78 0.36 0.04 65396.72\n",
" 9.75 3232.02 1366.12 0.30 0.04 62843.62\n",
" 10.05 2932.45 1239.49 0.23 0.03 59331.06\n",
" 10.35 2613.18 1104.54 0.19 0.02 54928.06\n",
" 10.65 2281.21 964.23 0.15 0.02 49741.36\n",
" 10.95 1944.28 821.81 0.11 0.01 43916.61\n",
" 11.25 1610.69 680.81 0.08 0.01 37638.04\n",
" 11.55 1289.19 544.92 0.06 0.01 31126.45\n",
" 11.85 988.66 417.89 0.04 0.00 24634.66\n",
" 12.15 717.91 303.45 0.02 0.00 18439.92\n",
" 12.45 485.17 205.07 0.01 0.00 12832.42\n",
" 12.75 297.73 125.84 0.01 0.00 8100.31\n",
" 13.05 161.08 68.09 0.00 0.00 4503.77\n",
" 13.35 75.82 32.05 0.00 0.00 2176.69\n",
" 13.65 31.84 13.46 0.00 0.00 937.66\n",
" 13.95 12.27 5.19 0.00 0.00 370.54\n",
" 14.25 4.31 1.82 0.00 0.00 133.28\n",
" 14.55 1.32 0.56 0.00 0.00 41.80\n",
"Totals 96108.76 40623.44 52.85 22.50 1284458.65\n"
]
}
],
"source": [
"# initialize counters\n",
"tot_nue_e = 0.\n",
"tot_nuebar_e = 0.\n",
"tot_nue_e_vis = 0.\n",
"tot_nuebar_e_vis = 0.\n",
"tot_nuebar_IBD = 0.\n",
"\n",
"print \"{:>14}{:>14}{:>14}{:>14}{:>14}{:>14}\". \\\n",
" format('enu(MeV)','n_nue_e','n_nuebar_e','n_nue_e_evis', \\\n",
" 'n_nuebar_e_evis','n_nuebar_IBD')\n",
"for ib in range (0,nebins-1):\n",
" de_ib = emax_ib/nebins\n",
" enu = (ib+0.5)*de_ib\n",
" line=\"{:14.2f}{:14.2f}{:14.2f}{:14.2f}{:14.2f}{:14.2f}\". \\\n",
" format(enu,nevts_nue_e[ib],nevts_nuebar_e[ib],nvis_nue_e[ib], \\\n",
" nvis_nuebar_e[ib],nevts_nuebar_IBD[ib])\n",
" print line\n",
" calcrunfile.write(line + \"\\n\")\n",
" tot_nue_e = tot_nue_e + nevts_nue_e[ib]\n",
" tot_nuebar_e = tot_nuebar_e + nevts_nuebar_e[ib]\n",
" tot_nuebar_IBD = tot_nuebar_IBD + nevts_nuebar_IBD[ib]\n",
" tot_nue_e_vis = tot_nue_e_vis + nvis_nue_e[ib]\n",
" tot_nuebar_e_vis = tot_nuebar_e_vis + nvis_nuebar_e[ib]\n",
" \n",
"line = \"{:14}{:14.2f}{:14.2f}{:14.2f}{:14.2f}{:14.2f}\". \\\n",
" format(\"Totals\",tot_nue_e,tot_nuebar_e,tot_nue_e_vis, \\\n",
" tot_nuebar_e_vis,tot_nuebar_IBD)\n",
"print line\n",
"calcrunfile.write(line + \"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 298,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Events vs levt\n",
" 6.131 2909\n",
" 6.391 8246\n",
" 6.651 1.295e+04\n",
" 6.911 1.709e+04\n",
" 7.171 2.072e+04\n",
" 7.431 2.391e+04\n",
" 7.69 2.669e+04\n",
" 7.95 2.91e+04\n",
" 8.21 3.119e+04\n",
" 8.47 3.297e+04\n",
" 8.73 3.449e+04\n",
" 8.99 3.575e+04\n",
" 9.25 3.679e+04\n",
" 9.51 3.762e+04\n",
" 9.77 3.825e+04\n",
" 10.03 3.871e+04\n",
" 10.29 3.9e+04\n",
" 10.55 3.914e+04\n",
" 10.81 3.914e+04\n",
" 11.07 3.901e+04\n",
" 11.33 3.875e+04\n",
" 11.59 3.839e+04\n",
" 11.85 3.791e+04\n",
" 12.11 3.734e+04\n",
" 12.37 3.667e+04\n",
" 12.63 3.592e+04\n",
" 12.89 3.508e+04\n",
" 13.15 3.416e+04\n",
" 13.41 3.318e+04\n",
" 13.67 3.212e+04\n",
" 13.93 3.1e+04\n",
" 14.19 2.982e+04\n",
" 14.45 2.858e+04\n",
" 14.71 2.728e+04\n",
" 14.97 2.594e+04\n",
" 15.23 2.454e+04\n",
" 15.49 2.31e+04\n",
" 15.75 2.161e+04\n",
" 16.01 2.008e+04\n",
" 16.27 1.851e+04\n",
" 16.53 1.69e+04\n",
" 16.79 1.526e+04\n",
" 17.05 1.357e+04\n",
" 17.31 1.186e+04\n",
" 17.57 1.011e+04\n",
" 17.83 8338\n",
" 18.09 6532\n",
" 18.35 4698\n",
" 18.61 2838\n",
" 18.87 713.2\n"
]
}
],
"source": [
"print \n",
"print 'Events vs levt'\n",
"calcrunfile.write(\"\\n\")\n",
"calcrunfile.write(\"Events vs levt\\n\")\n",
"for ib in range (0,nlbins):\n",
" dl_ib = (2.*detector_radius)/nlbins\n",
" levt = (ib+0.5)*dl_ib+(d_source_detector-detector_radius)\n",
" line = \"{:14.4g}{:14.4g}\".format(levt,nevts_nuebar_IBD_levt[ib])\n",
" print line\n",
" calcrunfile.write(line + \"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_In Python, any file that's been opened must be closed._"
]
},
{
"cell_type": "code",
"execution_count": 299,
"metadata": {},
"outputs": [],
"source": [
"calcrunfile.close()\n",
"eventfile.close()\n",
"evisfile.close()"
]
},
{
"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
}