{
"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": 1,
"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": 2,
"metadata": {},
"outputs": [],
"source": [
"nenu = 100 # Smimizu bins for integration\n",
"nl = 100 # distance bins for integration\n",
"ny = 100 # y bins for nuebar-e integration\n",
"\n",
"nebins = 50 # energy bins for printout\n",
"nlbins = 50 # distance bins for printout\n",
"emax_ib = 15.\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": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"10000 1000000 0.238\n"
]
}
],
"source": [
"eventfile = open(\"nuebar_event_array_Borexino.dat\",\"w\")\n",
"eventfile.write( str(nenu*nl) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"IBDfile = open(\"IBD_event_array_Borexino.dat\",\"w\")\n",
"IBDfile.write( str(nenu*nl)+ \" \" + \"enu,levt,xsecIBD,n_nuebar_IBD\" + \"\\n\")\n",
"\n",
"evisfile = open(\"nuebar_evis_array_Borexino.dat\",\"w\")\n",
"evisfile.write( str(nenu*nl*ny) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"nuebar_e_file = open(\"nuebar_e_event_array_Borexino.dat\",\"w\")\n",
"nuebar_e_file.write( str(nenu*nl*ny) + \" \" + str(sinsqthw) + \" \" + \"enu,levt,e_elect,xsec_nuebar_e,wtnuebar\" + \"\\n\")\n",
"\n",
"calcrunfile = open(\"current_istope_event_calc_run_Borexino.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": 4,
"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": 5,
"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: 100.0 radius: 3.01 shield: 5.0 distance: 14.00\n"
]
}
],
"source": [
"from math import pi\n",
"from random import random\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 = 100.0 # metric tons - use full volume of Borexino\n",
"\n",
"# detector_density = 0.78 # g/cm3 KamLAND\n",
"detector_density = 0.876 # g/cm3 Borexino\n",
"\n",
"electrons_per_molecule = 66.0 # Pseudocumine\n",
"mol_weight = 120.2 # Pseudocumine\n",
"free_protons_per_molecule = 12.0 # Pseudocumine\n",
"\n",
"detector_volume=(detector_mass*1000.)/(detector_density*1000.)\n",
"detector_radius = (.75*detector_mass*1000./(pi*detector_density*1000.))**(1./3.)\n",
"\n",
"d_source_detector = 14.0 # estimated distance from source to detector (m) Borexino\n",
"\n",
"lshield = d_source_detector-9.0 # 9m tank radius\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": 6,
"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": 7,
"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",
" #print float(e), float(flux)\n",
" \n",
"nshimi_file = len(tmp_e) # The number of entries in the shimizu file that pass the cut\n"
]
},
{
"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": 8,
"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": 9,
"metadata": {},
"outputs": [],
"source": [
"# rebin into nenu energy bins\n",
"nshimi = nshimi_file/nenu\n",
"#print 'nshimi_file,nenu,nshimi: ',nshimi_file,nenu,nshimi\n",
"#print 'Shimizu flux dist'\n",
"#for ienu in range (0,nenu-1):\n",
"for ienu in range (0,nenu): # range fuction stops before upper number\n",
" #for ib in range (ienu*nshimi, (ienu+1)*nshimi-1):\n",
" for ib in range (ienu*nshimi, (ienu+1)*nshimi): # range fuction stops before upper number\n",
" shimi_e[ienu] = shimi_e[ienu] + tmp_e[ib]*tmp_flux[ib]\n",
" shimi_flux[ienu] = shimi_flux[ienu] + tmp_flux[ib]\n",
" #print 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])\n",
" #print 'limits: ',ienu*nshimi, (ienu+1)*nshimi-1"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" enu(MeV) wt_spectrum flux0 IBD_xsec(enu)\n",
" 0.075 2.86912e-05 3.29182e+18 0\n",
" 0.225 0.000147516 1.6925e+19 0\n",
" 0.375 0.000372457 4.2733e+19 0\n",
" 0.525 0.000692596 7.94634e+19 0\n",
" 0.675 0.00109698 1.2586e+20 0\n",
" 0.825 0.00157584 1.808e+20 0\n",
" 0.975 0.00211987 2.43218e+20 0\n",
" 1.125 0.00272064 3.12146e+20 0\n",
" 1.275 0.00337008 3.86658e+20 0\n",
" 1.425 0.00406082 4.65909e+20 0\n",
" 1.575 0.00478578 5.49086e+20 0\n",
" 1.725 0.00553848 6.35446e+20 8.81831e-46\n",
" 1.875 0.00631268 7.24271e+20 1.91289e-44\n",
" 2.025 0.00710266 8.14908e+20 4.0625e-44\n",
" 2.175 0.00790292 9.06724e+20 6.54354e-44\n",
" 2.325 0.00870844 9.99144e+20 9.36229e-44\n",
" 2.475 0.00951439 1.09161e+21 1.25247e-43\n",
" 2.625 0.0103163 1.18362e+21 1.60366e-43\n",
" 2.775 0.0111101 1.27469e+21 1.99034e-43\n",
" 2.925 0.0118917 1.36437e+21 2.41302e-43\n",
" 3.075 0.0126576 1.45224e+21 2.87221e-43\n",
" 3.225 0.0134043 1.53792e+21 3.36838e-43\n",
" 3.375 0.0141288 1.62104e+21 3.90197e-43\n",
" 3.525 0.0148282 1.70128e+21 4.4734e-43\n",
" 3.675 0.0154998 1.77833e+21 5.08308e-43\n",
" 3.825 0.0161411 1.85191e+21 5.73138e-43\n",
" 3.975 0.01675 1.92177e+21 6.41866e-43\n",
" 4.125 0.0173244 1.98767e+21 7.14524e-43\n",
" 4.275 0.0178625 2.04941e+21 7.91143e-43\n",
" 4.425 0.0183627 2.1068e+21 8.71753e-43\n",
" 4.575 0.0188236 2.15969e+21 9.5638e-43\n",
" 4.725 0.019244 2.20792e+21 1.04505e-42\n",
" 4.875 0.0196228 2.25138e+21 1.13778e-42\n",
" 5.025 0.0199591 2.28997e+21 1.2346e-42\n",
" 5.175 0.0202523 2.3236e+21 1.33552e-42\n",
" 5.325 0.0205018 2.35222e+21 1.44057e-42\n",
" 5.475 0.0207071 2.37579e+21 1.54975e-42\n",
" 5.625 0.0208682 2.39427e+21 1.66308e-42\n",
" 5.775 0.0209849 2.40766e+21 1.78057e-42\n",
" 5.925 0.0210573 2.41597e+21 1.90223e-42\n",
" 6.075 0.0210857 2.41922e+21 2.02808e-42\n",
" 6.225 0.0210704 2.41747e+21 2.1581e-42\n",
" 6.375 0.0210119 2.41075e+21 2.29232e-42\n",
" 6.525 0.0209108 2.39916e+21 2.43073e-42\n",
" 6.675 0.020768 2.38277e+21 2.57334e-42\n",
" 6.825 0.0205843 2.36169e+21 2.72015e-42\n",
" 6.975 0.0203607 2.33604e+21 2.87115e-42\n",
" 7.125 0.0200983 2.30593e+21 3.02634e-42\n",
" 7.275 0.0197984 2.27153e+21 3.18572e-42\n",
" 7.425 0.0194624 2.23297e+21 3.34928e-42\n",
" 7.575 0.0190916 2.19043e+21 3.51703e-42\n",
" 7.725 0.0186877 2.14409e+21 3.68894e-42\n",
" 7.875 0.0182523 2.09413e+21 3.86501e-42\n",
" 8.025 0.0177871 2.04076e+21 4.04523e-42\n",
" 8.175 0.017294 1.98419e+21 4.22959e-42\n",
" 8.325 0.0167749 1.92463e+21 4.41807e-42\n",
" 8.475 0.0162319 1.86233e+21 4.61067e-42\n",
" 8.625 0.0156669 1.7975e+21 4.80736e-42\n",
" 8.775 0.0150822 1.73042e+21 5.00814e-42\n",
" 8.925 0.0144799 1.66132e+21 5.21298e-42\n",
" 9.075 0.0138624 1.59047e+21 5.42187e-42\n",
" 9.225 0.013232 1.51814e+21 5.63479e-42\n",
" 9.375 0.0125911 1.44461e+21 5.85172e-42\n",
" 9.525 0.0119421 1.37016e+21 6.07263e-42\n",
" 9.675 0.0112876 1.29506e+21 6.29752e-42\n",
" 9.825 0.01063 1.21961e+21 6.52635e-42\n",
" 9.975 0.00997196 1.14411e+21 6.7591e-42\n",
" 10.12 0.00931594 1.06884e+21 6.99575e-42\n",
" 10.28 0.0086646 9.94113e+20 7.23628e-42\n",
" 10.42 0.00802049 9.20213e+20 7.48066e-42\n",
" 10.57 0.00738624 8.47444e+20 7.72886e-42\n",
" 10.72 0.0067644 7.76099e+20 7.98086e-42\n",
" 10.88 0.00615755 7.06473e+20 8.23663e-42\n",
" 11.03 0.00556819 6.38854e+20 8.49614e-42\n",
" 11.17 0.0049988 5.73526e+20 8.75938e-42\n",
" 11.32 0.00445178 5.10765e+20 9.0263e-42\n",
" 11.47 0.00392947 4.50839e+20 9.29687e-42\n",
" 11.62 0.00343409 3.94003e+20 9.57108e-42\n",
" 11.78 0.00296776 3.40499e+20 9.84889e-42\n",
" 11.92 0.00253243 2.90553e+20 1.01303e-41\n",
" 12.07 0.00212992 2.44372e+20 1.04152e-41\n",
" 12.22 0.00176181 2.02138e+20 1.07036e-41\n",
" 12.38 0.0014295 1.64011e+20 1.09955e-41\n",
" 12.53 0.00113411 1.30119e+20 1.12908e-41\n",
" 12.67 0.000876499 1.00563e+20 1.15896e-41\n",
" 12.82 0.000657198 7.54021e+19 1.18917e-41\n",
" 12.97 0.000476267 5.46435e+19 1.21972e-41\n",
" 13.12 0.000332903 3.81949e+19 1.2506e-41\n",
" 13.28 0.000224653 2.5775e+19 1.28181e-41\n",
" 13.42 0.000147035 1.68698e+19 1.31335e-41\n",
" 13.57 9.38161e-05 1.07638e+19 1.3452e-41\n",
" 13.72 5.86341e-05 6.72726e+18 1.37738e-41\n",
" 13.88 3.59068e-05 4.11968e+18 1.40987e-41\n",
" 14.03 2.15348e-05 2.47075e+18 1.44267e-41\n",
" 14.17 1.25839e-05 1.44379e+18 1.47578e-41\n",
" 14.32 7.13389e-06 8.18491e+17 1.50919e-41\n",
" 14.47 3.88685e-06 4.45949e+17 1.54291e-41\n",
" 14.62 2.01884e-06 2.31627e+17 1.57692e-41\n",
" 14.77 9.83975e-07 1.12894e+17 1.61124e-41\n",
" 14.92 4.43241e-07 5.08543e+16 1.64584e-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",
"print \"{:>12}{:>14}{:<14}{:>14}\". \\\n",
" format('enu(MeV)','wt_spectrum',' flux0','IBD_xsec(enu)')\n",
"#for ienu in range (0,nenu-1):\n",
"for ienu in range (0,nenu): # range fuction stops before upper number\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)*free_protons_per_molecule # use parameterization with free protons per molecule\n",
" if enu <= 1.29332: \n",
" xsecIBD = 0.0\n",
"\n",
" print \"{:12.4g}{:14.6g}{:14.6g}{:14.6g}\".format(enu,wt_spectrum,flux0,IBD_xsec(enu))\n",
" \n",
" # nu-e scattering xsecs\n",
" eminus = -0.5-sinsqthw\n",
" eplus = -sinsqthw\n",
" # (use num of electrons and mol weight for detector material) Use grams and cm units.\n",
" # sig0 = 2Gmu^2*me*enu/(pi) in cm^2 (Should be 1.72328d-44 in Fortran before was 1.5e-44)\n",
" sig0 = 1.72328e-44*enu # enu in MeV\n",
" xsec_nue_e = sig0*(eminus**2+eplus**2/3.0)*electrons_per_molecule \n",
" xsec_nuebar_e = sig0*(eplus**2+eminus**2/3.0)*electrons_per_molecule\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",
" for il in range (1,nl+1): # range fuction stops before upper number\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",
" # (use num of electrons and mol weight for detector material) Use grams and cm units.\n",
" n_per_xsec = detector_density*(100.*dl)*(1./mol_weight)*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) # ib should be between 0 and nebins-1\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) # ib should be between 0 and nlbins-1\n",
" nevts_nuebar_IBD_levt[ib] = nevts_nuebar_IBD_levt[ib] + n_nuebar_IBD\n",
"\n",
" # write out event array with etrue\n",
" line = \"{:14.4g}{:18.8g}{:18.8g}{:18.8g}\".format(enu,levt,n_nuebar_IBD,n_nuebar)\n",
" eventfile.write(line + \"\\n\")\n",
" #if n_nuebar_IBD > 0:\n",
" line = \"{:14.4g}{:18.8g}{:18.8g}{:18.8g}\".format(enu,levt,xsecIBD,n_nuebar_IBD)\n",
" IBDfile.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+1): # range fuction stops before upper number\n",
" y = (float(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(e_elect/de_ib) # ib should be between 0 and nebins-1\n",
" if (ib < 0 or ib > nebins-1):\n",
" print 'ib problem - ienu,il,iy,ib,e_elect,de_ib: ',ienu,il,iy,ib,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",
" #if (ienu == 10 and il == 50):\n",
" #print iy,y,dy,e_elect,ib,de_ib,n_nuebar,(eplus**2+eminus**2*(1-y)**2),norm,wtnuebar\n",
" # write out event array with evis\n",
" line = \"{:12.4g}{:18.8g}{:18.8g}{:18.8g}{:18.8g}\".format(enu,levt,e_elect,wtnue,wtnuebar)\n",
" evisfile.write(line + \"\\n\")\n",
" line = \"{:12.4g}{:18.8g}{:18.8g}{:18.8g}{:18.8g}\".format(enu,levt,e_elect,xsec_nuebar_e,wtnuebar)\n",
" nuebar_e_file.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": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" enu/evis(MeV) n_nuebar_e_enu n_nuebar_e_evis n_nuebar_IBD_enu\n",
" 0.15 0.02 457.11 0.00\n",
" 0.45 0.32 412.13 0.00\n",
" 0.75 1.30 362.65 0.00\n",
" 1.05 3.27 329.16 0.00\n",
" 1.35 6.43 295.53 0.00\n",
" 1.65 10.90 261.20 0.14\n",
" 1.95 16.73 229.14 11.57\n",
" 2.25 23.89 208.15 37.66\n",
" 2.55 32.30 183.47 80.43\n",
" 2.85 41.86 165.30 143.59\n",
" 3.15 52.41 146.07 230.35\n",
" 3.45 63.77 129.87 343.27\n",
" 3.75 75.73 113.58 484.11\n",
" 4.05 88.08 100.27 653.68\n",
" 4.35 100.56 89.22 851.78\n",
" 4.65 112.96 79.82 1077.14\n",
" 4.95 125.03 68.41 1327.38\n",
" 5.25 136.52 63.35 1599.07\n",
" 5.55 147.23 51.33 1887.76\n",
" 5.85 156.92 46.67 2188.03\n",
" 6.15 165.41 41.30 2493.66\n",
" 6.45 172.52 35.95 2797.73\n",
" 6.75 178.08 30.38 3092.80\n",
" 7.05 181.98 26.66 3371.09\n",
" 7.35 184.10 23.44 3624.72\n",
" 7.65 184.38 20.12 3845.90\n",
" 7.95 182.78 16.75 4027.19\n",
" 8.25 179.30 13.99 4161.75\n",
" 8.55 173.98 12.06 4243.62\n",
" 8.85 166.89 10.05 4267.95\n",
" 9.15 158.15 8.26 4231.29\n",
" 9.45 147.89 6.75 4131.81\n",
" 9.75 136.31 5.03 3969.57\n",
" 10.05 123.65 4.23 3746.71\n",
" 10.35 110.15 3.32 3467.61\n",
" 10.65 96.12 2.34 3139.07\n",
" 10.95 81.89 1.84 2770.34\n",
" 11.25 67.81 1.26 2373.09\n",
" 11.55 54.24 0.87 1961.33\n",
" 11.85 41.56 0.57 1551.08\n",
" 12.15 30.15 0.35 1159.88\n",
" 12.45 20.35 0.19 806.10\n",
" 12.75 12.47 0.10 507.96\n",
" 13.05 6.73 0.05 281.83\n",
" 13.35 3.16 0.02 135.96\n",
" 13.65 1.33 0.01 58.49\n",
" 13.95 0.51 0.00 23.09\n",
" 14.25 0.18 0.00 8.29\n",
" 14.55 0.05 0.00 2.59\n",
" 14.85 0.01 0.00 0.65\n",
"Totals 4058.40 4058.32 81169.10\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",
"#line = \"{:>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",
"\n",
"line = \"{:>18}{:>18}{:>18}{:>18}\". \\\n",
" format('enu/evis(MeV)','n_nuebar_e_enu', \\\n",
" 'n_nuebar_e_evis','n_nuebar_IBD_enu')\n",
"\n",
"print line\n",
"calcrunfile.write(line + \"\\n\")\n",
"\n",
"for ib in range (0,nebins): # range fuction stops before upper number\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",
" line=\"{:18.2f}{:18.2f}{:18.2f}{:18.2f}\". \\\n",
" format(enu,nevts_nuebar_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 = \"{:18}{:18.2f}{:18.2f}{:18.2f}\". \\\n",
" format(\"Totals\",tot_nuebar_e, \\\n",
" tot_nuebar_e_vis,tot_nuebar_IBD)\n",
"print line\n",
"calcrunfile.write(line + \"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" L_event(m) IBD Events\n",
" 11.05 120.55\n",
" 11.17 351.43\n",
" 11.29 567.81\n",
" 11.41 770.16\n",
" 11.53 958.92\n",
" 11.65 1134.50\n",
" 11.77 1297.31\n",
" 11.89 1447.74\n",
" 12.01 1586.16\n",
" 12.13 1712.93\n",
" 12.25 1828.39\n",
" 12.37 1932.87\n",
" 12.50 2026.68\n",
" 12.62 2110.14\n",
" 12.74 2183.54\n",
" 12.86 2247.16\n",
" 12.98 2301.27\n",
" 13.10 2346.13\n",
" 13.22 2382.00\n",
" 13.34 2409.12\n",
" 13.46 2427.73\n",
" 13.58 2438.05\n",
" 13.70 2440.30\n",
" 13.82 2434.68\n",
" 13.94 2421.42\n",
" 14.06 2400.69\n",
" 14.18 2372.69\n",
" 14.30 2337.61\n",
" 14.42 2295.62\n",
" 14.54 2246.90\n",
" 14.66 2191.60\n",
" 14.78 2129.90\n",
" 14.90 2061.94\n",
" 15.02 1987.87\n",
" 15.14 1907.85\n",
" 15.26 1822.01\n",
" 15.38 1730.48\n",
" 15.50 1633.41\n",
" 15.63 1530.92\n",
" 15.75 1423.13\n",
" 15.87 1310.17\n",
" 15.99 1192.14\n",
" 16.11 1069.17\n",
" 16.23 941.37\n",
" 16.35 808.84\n",
" 16.47 671.68\n",
" 16.59 530.00\n",
" 16.71 383.89\n",
" 16.83 233.45\n",
" 16.95 78.77\n"
]
}
],
"source": [
"print \n",
"print ' L_event(m) IBD Events'\n",
"calcrunfile.write(\"\\n\")\n",
"calcrunfile.write(\" L_event(m) IBD Events\\n\")\n",
"for ib in range (0,nlbins): # range fuction stops before upper number\n",
" dl_ib = (2.*detector_radius)/nlbins\n",
" levt = (ib+0.5)*dl_ib+(d_source_detector-detector_radius)\n",
" line = \"{:14.2f}{:16.2f}\".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": 13,
"metadata": {},
"outputs": [],
"source": [
"calcrunfile.close()\n",
"eventfile.close()\n",
"evisfile.close()\n",
"nuebar_e_file.close()\n",
"IBDfile.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
}