{
"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_SuperK.dat\",\"w\")\n",
"eventfile.write( str(nenu*nl) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"IBDfile = open(\"IBD_event_array_SuperK.dat\",\"w\")\n",
"IBDfile.write( str(nenu*nl)+ \" \" + \"enu,levt,xsecIBD,n_nuebar_IBD\" + \"\\n\")\n",
"\n",
"evisfile = open(\"nuebar_evis_array_SuperK.dat\",\"w\")\n",
"evisfile.write( str(nenu*nl*ny) + \" \" + str(sinsqthw) + \"\\n\")\n",
"\n",
"nuebar_e_file = open(\"nuebar_e_event_array_SuperK.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_SuperK.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",
"duty_fact = 0.8\n",
"nyrs = 5.0*duty_fact # 5yrs with DF\n",
"\n",
"eff_det = 0.69 # pure water\n",
"eff_Gd_tag = 0.69 # pure water neutron capture\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 (MeV): 15.000 eff_det: 0.690 eff_Gd_tag: 0.690\n",
"mass: 22500.0 height: 32.20 radius: 14.91 shield: 9.25 distance: 9.25\n"
]
}
],
"source": [
"from math import pi\n",
"from random import random\n",
"from math import sqrt # In python, this is not built-in.\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",
"# 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",
"detector_mass = 22500. # metric tons - SuperK with 2m cut from ID edges\n",
"\n",
"# detector_density = 0.85 # g/cm3\n",
"detector_density = 1.00 # g/cm3\n",
"\n",
"electrons_per_molecule = 10.0 # water\n",
"mol_weight = 18. # water\n",
"free_protons_per_molecule = 2.0 # water\n",
"\n",
"# inner detector: height = 36.2m - 4m = 32.2m; radius = 16.9m - 2m = 14.9m\n",
"# below 4.5MeV there is a tighter fid_cut (see 1606.07538)\n",
"detector_height = 32.2 # SuperK height of cylinder in m\n",
"\n",
"detector_volume=(detector_mass*1000.)/(detector_density*1000.)\n",
"detector_radius = sqrt(detector_volume/(pi*detector_height))\n",
"\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",
"lshield = 9.25 # IsoDAR_shield+water_buffer+fid_cut=4.5m+2.75m+2m\n",
"d_source_detector = lshield\n",
"\n",
"line = \"beta-decay endpoint (MeV): {:6.3f} eff_det: {:5.3f} eff_Gd_tag: {:5.3f}\".format(e_endpoint,eff_det,eff_Gd_tag)\n",
"calcrunfile.write(line + \"\\n\")\n",
"print line\n",
"\n",
"line = \"mass: {:8.1f} height: {:6.2f} radius: {:6.2f} shield: {:6.2f} distance: {:6.2f}\".format(detector_mass,detector_height,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 2.27136e+18 0\n",
" 0.225 0.000147516 1.16782e+19 0\n",
" 0.375 0.000372457 2.94858e+19 0\n",
" 0.525 0.000692596 5.48298e+19 0\n",
" 0.675 0.00109698 8.68435e+19 0\n",
" 0.825 0.00157584 1.24752e+20 0\n",
" 0.975 0.00211987 1.67821e+20 0\n",
" 1.125 0.00272064 2.15381e+20 0\n",
" 1.275 0.00337008 2.66794e+20 0\n",
" 1.425 0.00406082 3.21477e+20 0\n",
" 1.575 0.00478578 3.7887e+20 0\n",
" 1.725 0.00553848 4.38457e+20 8.81831e-46\n",
" 1.875 0.00631268 4.99747e+20 1.91289e-44\n",
" 2.025 0.00710266 5.62286e+20 4.0625e-44\n",
" 2.175 0.00790292 6.2564e+20 6.54354e-44\n",
" 2.325 0.00870844 6.89409e+20 9.36229e-44\n",
" 2.475 0.00951439 7.53212e+20 1.25247e-43\n",
" 2.625 0.0103163 8.16699e+20 1.60366e-43\n",
" 2.775 0.0111101 8.79535e+20 1.99034e-43\n",
" 2.925 0.0118917 9.41413e+20 2.41302e-43\n",
" 3.075 0.0126576 1.00204e+21 2.87221e-43\n",
" 3.225 0.0134043 1.06116e+21 3.36838e-43\n",
" 3.375 0.0141288 1.11852e+21 3.90197e-43\n",
" 3.525 0.0148282 1.17388e+21 4.4734e-43\n",
" 3.675 0.0154998 1.22705e+21 5.08308e-43\n",
" 3.825 0.0161411 1.27782e+21 5.73138e-43\n",
" 3.975 0.01675 1.32602e+21 6.41866e-43\n",
" 4.125 0.0173244 1.37149e+21 7.14524e-43\n",
" 4.275 0.0178625 1.41409e+21 7.91143e-43\n",
" 4.425 0.0183627 1.4537e+21 8.71753e-43\n",
" 4.575 0.0188236 1.49018e+21 9.5638e-43\n",
" 4.725 0.019244 1.52346e+21 1.04505e-42\n",
" 4.875 0.0196228 1.55345e+21 1.13778e-42\n",
" 5.025 0.0199591 1.58008e+21 1.2346e-42\n",
" 5.175 0.0202523 1.60329e+21 1.33552e-42\n",
" 5.325 0.0205018 1.62303e+21 1.44057e-42\n",
" 5.475 0.0207071 1.63929e+21 1.54975e-42\n",
" 5.625 0.0208682 1.65205e+21 1.66308e-42\n",
" 5.775 0.0209849 1.66128e+21 1.78057e-42\n",
" 5.925 0.0210573 1.66702e+21 1.90223e-42\n",
" 6.075 0.0210857 1.66926e+21 2.02808e-42\n",
" 6.225 0.0210704 1.66805e+21 2.1581e-42\n",
" 6.375 0.0210119 1.66342e+21 2.29232e-42\n",
" 6.525 0.0209108 1.65542e+21 2.43073e-42\n",
" 6.675 0.020768 1.64411e+21 2.57334e-42\n",
" 6.825 0.0205843 1.62957e+21 2.72015e-42\n",
" 6.975 0.0203607 1.61187e+21 2.87115e-42\n",
" 7.125 0.0200983 1.59109e+21 3.02634e-42\n",
" 7.275 0.0197984 1.56735e+21 3.18572e-42\n",
" 7.425 0.0194624 1.54075e+21 3.34928e-42\n",
" 7.575 0.0190916 1.5114e+21 3.51703e-42\n",
" 7.725 0.0186877 1.47942e+21 3.68894e-42\n",
" 7.875 0.0182523 1.44495e+21 3.86501e-42\n",
" 8.025 0.0177871 1.40813e+21 4.04523e-42\n",
" 8.175 0.017294 1.36909e+21 4.22959e-42\n",
" 8.325 0.0167749 1.328e+21 4.41807e-42\n",
" 8.475 0.0162319 1.28501e+21 4.61067e-42\n",
" 8.625 0.0156669 1.24028e+21 4.80736e-42\n",
" 8.775 0.0150822 1.19399e+21 5.00814e-42\n",
" 8.925 0.0144799 1.14631e+21 5.21298e-42\n",
" 9.075 0.0138624 1.09743e+21 5.42187e-42\n",
" 9.225 0.013232 1.04752e+21 5.63479e-42\n",
" 9.375 0.0125911 9.96783e+20 5.85172e-42\n",
" 9.525 0.0119421 9.45407e+20 6.07263e-42\n",
" 9.675 0.0112876 8.93592e+20 6.29752e-42\n",
" 9.825 0.01063 8.41533e+20 6.52635e-42\n",
" 9.975 0.00997196 7.89436e+20 6.7591e-42\n",
" 10.12 0.00931594 7.37502e+20 6.99575e-42\n",
" 10.28 0.0086646 6.85938e+20 7.23628e-42\n",
" 10.42 0.00802049 6.34947e+20 7.48066e-42\n",
" 10.57 0.00738624 5.84737e+20 7.72886e-42\n",
" 10.72 0.0067644 5.35508e+20 7.98086e-42\n",
" 10.88 0.00615755 4.87467e+20 8.23663e-42\n",
" 11.03 0.00556819 4.40809e+20 8.49614e-42\n",
" 11.17 0.0049988 3.95733e+20 8.75938e-42\n",
" 11.32 0.00445178 3.52428e+20 9.0263e-42\n",
" 11.47 0.00392947 3.11079e+20 9.29687e-42\n",
" 11.62 0.00343409 2.71862e+20 9.57108e-42\n",
" 11.78 0.00296776 2.34945e+20 9.84889e-42\n",
" 11.92 0.00253243 2.00481e+20 1.01303e-41\n",
" 12.07 0.00212992 1.68617e+20 1.04152e-41\n",
" 12.22 0.00176181 1.39475e+20 1.07036e-41\n",
" 12.38 0.0014295 1.13167e+20 1.09955e-41\n",
" 12.53 0.00113411 8.97824e+19 1.12908e-41\n",
" 12.67 0.000876499 6.93886e+19 1.15896e-41\n",
" 12.82 0.000657198 5.20275e+19 1.18917e-41\n",
" 12.97 0.000476267 3.7704e+19 1.21972e-41\n",
" 13.12 0.000332903 2.63545e+19 1.2506e-41\n",
" 13.28 0.000224653 1.77848e+19 1.28181e-41\n",
" 13.42 0.000147035 1.16401e+19 1.31335e-41\n",
" 13.57 9.38161e-05 7.42701e+18 1.3452e-41\n",
" 13.72 5.86341e-05 4.64181e+18 1.37738e-41\n",
" 13.88 3.59068e-05 2.84258e+18 1.40987e-41\n",
" 14.03 2.15348e-05 1.70482e+18 1.44267e-41\n",
" 14.17 1.25839e-05 9.96215e+17 1.47578e-41\n",
" 14.32 7.13389e-06 5.64759e+17 1.50919e-41\n",
" 14.47 3.88685e-06 3.07705e+17 1.54291e-41\n",
" 14.62 2.01884e-06 1.59823e+17 1.57692e-41\n",
" 14.77 9.83975e-07 7.7897e+16 1.61124e-41\n",
" 14.92 4.43241e-07 3.50894e+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*eff_det \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 correct 2 free protons per CH2\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",
" # 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 # electrons per molecule\n",
" xsec_nuebar_e = sig0*(eplus**2+eminus**2/3.0)*electrons_per_molecule # 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",
" ytop = detector_height+d_source_detector\n",
" lmax = sqrt(detector_radius**2+ytop**2)\n",
" dl = (lmax-d_source_detector)/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",
" levt = dl*(il-0.5)+d_source_detector\n",
"# region 1 where xevt is less than detector radius\n",
" yevt = d_source_detector\n",
" xevt = sqrt(levt**2-yevt**2)\n",
" if xevt <= detector_radius:\n",
"# costh = sqrt(1.-(yevt/levt)**2)\n",
" costh = sqrt(1.-(xevt/levt)**2) # change in geometry for cylinder\n",
"# region 2 where xevt equals detector radius\n",
" else:\n",
" xevt = detector_radius\n",
" yevt = sqrt(levt**2-xevt**2)\n",
"# costh = sqrt(1.-(yevt/levt)**2)\n",
" costh = sqrt(1.-(xevt/levt)**2) # change in geometry for cylinder\n",
"\n",
" # calculate the flux into this solid angle\n",
" flux = flux0*0.5*(1-costh)\n",
" \n",
"# need to subtract as solid angle sphere grows bigger than detector height\n",
" if levt > ytop:\n",
" ysub = ytop\n",
" xsub = sqrt(levt**2-ysub**2)\n",
"# costhsub = sqrt(1.-(ysub/levt)**2)\n",
" costhsub = sqrt(1.-(xsub/levt)**2) # change in geometry for cylinder\n",
" flux = flux - flux0*0.5*(1-costhsub)\n",
" if flux < 0.:\n",
" print levt,flux,xevt,yevt,costh,flux0*0.5*(1-costh)\n",
" print levt,flux,xsub,ysub,costhsub,flux0*0.5*(1-costhsub) \n",
" \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./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*eff_Gd_tag\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 = (lmax-d_source_detector)/nlbins\n",
" ib = int((levt-d_source_detector)/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\")\n",
" #end iy loop\n",
" #end il loop\n",
"#end ienu loop"
]
},
{
"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 1.27 25678.39 0.00\n",
" 0.45 18.04 23151.58 0.00\n",
" 0.75 73.14 20371.78 0.00\n",
" 1.05 183.78 18490.71 0.00\n",
" 1.35 361.41 16601.55 0.00\n",
" 1.65 612.58 14672.65 5.89\n",
" 1.95 939.73 12871.91 493.20\n",
" 2.25 1341.75 11692.89 1605.56\n",
" 2.55 1814.59 10306.58 3429.41\n",
" 2.85 2351.67 9285.70 6122.20\n",
" 3.15 2944.39 8205.72 9821.29\n",
" 3.45 3582.48 7295.54 14635.99\n",
" 3.75 4254.40 6380.20 20640.94\n",
" 4.05 4947.67 5632.88 27870.97\n",
" 4.35 5649.18 5012.12 36317.41\n",
" 4.65 6345.54 4483.94 45925.84\n",
" 4.95 7023.29 3842.92 56595.42\n",
" 5.25 7669.21 3558.82 68179.57\n",
" 5.55 8270.55 2883.65 80488.10\n",
" 5.85 8815.24 2621.91 93290.72\n",
" 6.15 9292.16 2320.04 106321.83\n",
" 6.45 9691.25 2019.38 119286.41\n",
" 6.75 10003.79 1706.82 131867.27\n",
" 7.05 10222.49 1497.35 143732.98\n",
" 7.35 10341.68 1316.65 154546.98\n",
" 7.65 10357.44 1130.02 163977.28\n",
" 7.95 10267.71 940.76 171706.71\n",
" 8.25 10072.44 785.66 177444.03\n",
" 8.55 9773.60 677.27 180934.62\n",
" 8.85 9375.31 564.72 181972.03\n",
" 9.15 8883.84 463.90 180408.85\n",
" 9.45 8307.65 379.27 176167.50\n",
" 9.75 7657.39 282.52 169250.18\n",
" 10.05 6945.82 237.53 159747.93\n",
" 10.35 6187.71 186.41 147848.18\n",
" 10.65 5399.75 131.64 133840.31\n",
" 10.95 4600.31 103.56 118118.56\n",
" 11.25 3809.13 70.53 101181.22\n",
" 11.55 3046.93 48.85 83625.08\n",
" 11.85 2334.86 31.88 66133.13\n",
" 12.15 1693.75 19.80 49453.76\n",
" 12.45 1143.15 10.60 34369.74\n",
" 12.75 700.27 5.64 21657.67\n",
" 13.05 378.09 2.70 12016.56\n",
" 13.35 177.64 1.16 5796.80\n",
" 13.65 74.49 0.45 2493.86\n",
" 13.95 28.68 0.16 984.36\n",
" 14.25 10.06 0.05 353.51\n",
" 14.55 3.07 0.01 110.62\n",
" 14.85 0.76 0.00 27.89\n",
"Totals 227981.10 227976.76 3460798.33\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",
" 8.90 13126.14\n",
" 9.60 36945.75\n",
" 10.29 57745.83\n",
" 10.99 76066.43\n",
" 11.69 92326.04\n",
" 12.38 106853.99\n",
" 13.08 119912.97\n",
" 13.77 131715.08\n",
" 14.47 142433.47\n",
" 15.17 152210.86\n",
" 15.86 161165.99\n",
" 16.56 169398.40\n",
" 17.25 162072.71\n",
" 17.95 146075.99\n",
" 18.65 132708.96\n",
" 19.34 121335.94\n",
" 20.04 111526.55\n",
" 20.73 102973.91\n",
" 21.43 95450.66\n",
" 22.13 88783.47\n",
" 22.82 82837.06\n",
" 23.52 77503.90\n",
" 24.21 72697.20\n",
" 24.91 68345.92\n",
" 25.61 64391.33\n",
" 26.30 60784.34\n",
" 27.00 57483.63\n",
" 27.69 54454.07\n",
" 28.39 51665.67\n",
" 29.09 49092.60\n",
" 29.78 46712.55\n",
" 30.48 44506.10\n",
" 31.17 42456.31\n",
" 31.87 40548.28\n",
" 32.57 38768.93\n",
" 33.26 37106.66\n",
" 33.96 35551.20\n",
" 34.66 34093.40\n",
" 35.35 32725.09\n",
" 36.05 31438.95\n",
" 36.74 30228.42\n",
" 37.44 29087.59\n",
" 38.14 28011.11\n",
" 38.83 26994.16\n",
" 39.53 26032.37\n",
" 40.22 25121.76\n",
" 40.92 22776.17\n",
" 41.62 16011.75\n",
" 42.31 9443.84\n",
" 43.01 3098.86\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 = (lmax-d_source_detector)/nlbins\n",
" levt = (ib-0.5)*dl_ib+d_source_detector\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
}