In [30]:
#Imports
from math import *
import numpy as np
import scipy as sp
import scipy.special
import scipy.interpolate as interpolate
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cbook as cbook
import seaborn as sns
import sys
import os
#Import custom modules
from physics import *
%matplotlib notebook
In [31]:
B0 = 1710.0/10**4 #Magnetic field strength in Tesla
def KEcalc(z,y):
"""Returns KE in J given z-position in m"""
return me*c**2*(sqrt((q*B0/(me*c))**2*((z**2+y**2)/(2*y))**2+1)-1)
def Radius(KE):
"""#Radius of electron orbit in m given KE in keV"""
return me*c/(q*B0)*sqrt((KE*1000*q/(me*c**2)+1)**2-1)
def zfcalc(KE,y):
"""#Returns z-position at screen in inches given KE in keV"""
R = Radius(KE)
return sqrt(R**2 - (y-R)**2)
def zfcalcGeneral(KE,yM,y):
R = Radius(KE)
zM = zfcalc(KE,yM)
return zM + (y - yM)*(R - yM)/zM
def KEcalcGeneral(zf,yM,yf):
"""Returns KE in J given z-position of electrons, y-position of magnet edge, and y-position of screen, all in m"""
a = (yM+yf)**2
b = -2*yM*(yf*(yM+yf)+zf**2)
d = yM**2*(zf**2+yf**2)
f = (me*c)/(q*B0)
g = (-b+sqrt(b**2-4*a*d))/(2*a)
return me*c**2*(sqrt(g**2+f**2)/f - 1)
def AngleIncidence(KE,yM):
R = Radius(KE)
return asin((R-yM)/R)
In [32]:
def getfns(folder,ext=''):
"""Get a list of full path filenames for all files in a folder and subfolders for the given extension"""
fns = []
for file in os.listdir(folder):
if file.endswith(ext):
fns.append(os.path.join(folder,file))
return fns
def readcsv(filename):
"""Read in a csv file and load it as an array"""
return np.loadtxt(open(filename, "rb"), delimiter=",")
def readcsvs(filenames):
"""Read in multiple csv files and load them in to a two-dimensional array"""
template = readcsv(filenames[0])
numfiles = len(filenames)
Data = np.zeros([numfiles,len(template)])
for i in range(numfiles):
spectrum = readcsv(filenames[i])
Data[i,:] = spectrum
return Data
def DataClean(Data):
"""Read in data and clean it by removing rows that have saturated pixel values"""
maxes = np.max(Data[:,500:],1)
includes = maxes<(2**16-1)
rejects = maxes>(2**16-2)
CleanData = Data[includes,:]
return CleanData
def DataAverage(Data):
"""Average input 2d array into a 1d array"""
return np.mean(Data,0)
In [33]:
MagnetOutFolderPath = os.curdir + '/Data/2015-08-17_Tiger_stripes_test/even_more_no_magnet'
MagnetInFolderPath = os.curdir + '/Data/2015-08-17_Tiger_stripes_test/magnet_in'
NoLaserFolderPath = os.curdir + '/Data/2015-08-17_Tiger_stripes_test/no_laser'
MagnetOutFiles = getfns(MagnetOutFolderPath,'')
MagnetInFiles = getfns(MagnetInFolderPath,'')
NoLaserFiles = getfns(NoLaserFolderPath,'')
MagnetOutData = readcsvs(MagnetOutFiles)
MagnetInData = readcsvs(MagnetInFiles)
Background = DataAverage(MagnetInData)
MagnetOutDataClean = DataClean(MagnetOutData)-Background
In [34]:
sns.set(font_scale=1.5)
fig1 = plt.figure(figsize=(8,6))
ax1 = fig1.add_subplot(111)
ax1.set_xlim(500,3648)
ax1.set_ylim(0,len(MagnetInData))
ax1.set_xlabel('Pixel')
ax1.set_ylabel('Instance')
# mesh1 = ax1.pcolormesh(MagnetInData, cmap='hot',vmin=0, vmax=10000)
mesh1 = ax1.pcolormesh(MagnetOutDataClean, cmap='inferno',vmin=0, vmax=65535)
In [35]:
MagnetOutDataCleanAvg = DataAverage(MagnetOutDataClean)
MagnetOutDataCleanAvgHat = savitzky_golay(MagnetOutDataCleanAvg,51,3) #Smoothed Data
fig2 = plt.figure(figsize=(12,8))
ax2 = fig2.add_subplot(111)
ax2.set_xlim(0,3648)
ax2.set_ylim(0,65535)
ax2.set_xlabel('Pixel Number')
ax2.set_ylabel('Pixel Value')
#ax2.semilogy()
ax2.plot(MagnetOutDataCleanAvg,linewidth=1)
ax2.plot(MagnetOutDataCleanAvgHat,linewidth=1,color='r')
Out[35]:
In [36]:
#Pixel values were determined by picking a region within each peak and trough
LanexPixelArray = np.hstack((np.arange(240,260),np.arange(670,750),np.arange(1090,1300),np.arange(1600,1700),\
np.arange(1930,2030),np.arange(2400,2450),np.arange(2820,2920),np.arange(3280,3350)))
LanexData = np.hstack((MagnetOutDataCleanAvgHat[240:260],MagnetOutDataCleanAvgHat[670:750],\
MagnetOutDataCleanAvgHat[1090:1300],MagnetOutDataCleanAvgHat[1600:1700],\
MagnetOutDataCleanAvgHat[1930:2030],MagnetOutDataCleanAvgHat[2400:2450],\
MagnetOutDataCleanAvgHat[2820:2920],MagnetOutDataCleanAvgHat[3280:3350]))
VinylPixelArray = np.hstack((np.arange(460,490),np.arange(900,950),np.arange(1400,1500),np.arange(1750,1850),\
np.arange(2150,2250),np.arange(2540,2640)))
VinylData = np.hstack((MagnetOutDataCleanAvgHat[460:490],MagnetOutDataCleanAvgHat[900:950],\
MagnetOutDataCleanAvgHat[1400:1500],MagnetOutDataCleanAvgHat[1750:1850],\
MagnetOutDataCleanAvgHat[2150:2250],MagnetOutDataCleanAvgHat[2540:2640]))
LanexInterpFunc = interpolate.interp1d(LanexPixelArray,LanexData,kind='slinear')
LanexPixelArrayInterp = np.arange(240,3350)
LanexDataInterp = LanexInterpFunc(LanexPixelArrayInterp)
VinylInterpFunc = interpolate.interp1d(VinylPixelArray,VinylData,kind='slinear')
VinylPixelArrayInterp = np.arange(460,2640)
VinylDataInterp = VinylInterpFunc(VinylPixelArrayInterp)
In [37]:
# sns.set(context='poster',font_scale=1.5)
# sns.set_style("darkgrid")
# sns.set_palette(palette='deep')
# sns.set_color_codes(palette='deep')
plt.figure(figsize=(12,6))
#With Lanex
#plt.plot(LanexPixelArray,LanexData,linestyle='None',marker='.')
plt.plot(LanexPixelArrayInterp,LanexDataInterp, linewidth=2,linestyle='--',label='Without Vinyl',color='b')
#Without Lanex
#plt.plot(VinylPixelArray,VinylData,linestyle='None',marker='.')
plt.plot(VinylPixelArrayInterp,VinylDataInterp,linewidth=2,linestyle='--',label='With Vinyl',color='r')
#All smoothed data
plt.plot(MagnetOutDataCleanAvgHat,linewidth=2,label='All Data',color='g')
plt.xlim(0,2000)
plt.xlabel('Pixel Number')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.ylabel('Pixel Value')
plt.legend()
plt.subplots_adjust(left=0.12,bottom=0.14) #Adjust spacing to prevent clipping of x and y labels
plt.savefig('strip_test.svg')
In [38]:
yM = 0.5 #magnet edge position in inches
CCDpos = 3.02 #CCD y-position relative to magnet edge in mm
yM = yM*.0254 #Convert to base units
CCDpos = CCDpos*10**-3 #Convert to base units
yf = yM + CCDpos #Screen position
KEcalcGeneralVec = np.vectorize(KEcalcGeneral) #Vectorize espec function for KE
ChosenPixels = [479.37419602, 658.07909036, 805.53608568, 929.62974462, 1079.52432829, 1325.661891,\
1454.37508809, 1637.97442378, 1801.34325342, 1984.12923212, 2202.34572023, 2412.63089924,\
2595.3896852]
ChosenPositions = np.multiply(np.add(np.multiply(ChosenPixels,8.0/10**3),9.5),10**-3) #Convert to meters
ChosenEnergies = KEcalcGeneralVec(ChosenPositions,yM,yf) #Energies corresponding to the chosen pixels
SignalRatio = LanexInterpFunc(ChosenPixels)/VinylInterpFunc(ChosenPixels)
print(ChosenEnergies/(1000*q),'\n')
print(SignalRatio)
MCNP was ran at the above chosen energies both with and without vinyl covering the lanex. The CCD depletion region thickness was also varied from 1 um to 20 um for each energy. The signal ratio at each of these energies can be determined by finding the ratio between the total signal without vinyl and the total signal with vinyl. This can then be compared with the experimentally determined signal ratio in order to find the appropriate depletion region thickness.
First, we need to read in the MCNP data for the "without vinyl" case. The energy deposited into each segment of the phosphor layer can be added to an array (code adapted from "Read MCNP Output" script):
</font>
In [39]:
mass = 2.95680*10**-2 #mass of active layer in g
Nsegments = 101 #number of segments
segmentmass = mass/Nsegments #mass of each segment in g
phosphorthickness = 81.*10**-6 #thickness of active layer in m
KEsim = np.multiply([338, 391, 440, 485],10**-3) #Chosen energies (only the first 4 were used to save time)
#Create array of angles
CosThetaSim = []
for x in KEsim:
CosThetaSim.append(np.cos(AngleIncidence(x*1000,yM)))
CCDthickness = np.arange(1,21,1) #Depletion region thickness array
MCNP_Directory = '/home/drake/Documents/Physics/Research/Python/MCNP_Code'
#Read files. As far as visible photons are concerned, the depletion region thickness doesn't matter
deptharray= []
Eabsorbed_Segments = []
Eabsorbed_Segments_Error = []
for i in range(len(KEsim)):
deptharray.append([])
Eabsorbed_Segments.append([])
Eabsorbed_Segments_Error.append([])
#Directory of the input files where vinyl was not used in the simulation
if i%2==0: #Energies where vinyl was experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output2/Out_{KE}MeV_{Theta}Degrees_20umCCD_inverted'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)))
else: #Energies where vinyl wasn't experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output1/Out_{KE}MeV_{Theta}Degrees_20umCCD'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)))
for segmentnumber in range(Nsegments-1):
segment = segmentnumber+100 #segment number
deptharray[i].append(segmentnumber*phosphorthickness/Nsegments) #depth for each segment in microns
printflag = False
with open(directory) as searchfile:
for line in searchfile:
left,sep,right = line.partition(' -{segmentlabel} '.format(segmentlabel=str(segment)))
if printflag:
Eabs_per_g = float(line[17:28])
Eabs_per_g_Error = float(line[29:35])
Eabsorbed_Segments[i].append(Eabs_per_g*segmentmass)
Eabsorbed_Segments_Error[i].append(Eabs_per_g_Error*segmentmass)
printflag = False
if sep:
printflag = True
Plot energy absorbed versus segment depth to make sure everything looks good:
In [40]:
plt.figure(figsize=(12,6))
plt.plot(np.multiply(deptharray[1],10**6),np.multiply(Eabsorbed_Segments[0],10**3),label='338 keV')
plt.plot(np.multiply(deptharray[1],10**6),np.multiply(Eabsorbed_Segments[1],10**3),label='391 keV')
plt.plot(np.multiply(deptharray[1],10**6),np.multiply(Eabsorbed_Segments[2],10**3),label='440 keV')
plt.plot(np.multiply(deptharray[1],10**6),np.multiply(Eabsorbed_Segments[3],10**3),label='485 keV')
plt.xlabel('Depth (um)')
plt.ylabel('Energy absorbed (keV)')
plt.legend()
Out[40]:
Next, the number of photons that reach the CCD can be calculated. From the quantum efficiency of the CCD, the contribution to the signal from visible photons can be found.
In [41]:
scatteringlength = 2.84*10**-6 #photon scattering length in Gd2O2S in m
conversionefficiency = 0.16 #electron energy to light energy conversion efficiency
emissionwavelength = 545*10**-9 #lanex emission wavelength in m
emissionenergy = hbar*2*pi*c/emissionwavelength
photonNumberArray = []
for i in range(len(KEsim)):
photonNumber = 0
for j in range(len(deptharray[i])):
photonNumber += conversionefficiency*Eabsorbed_Segments[i][j]*10**6*q/emissionenergy\
*(j+0.5)/Nsegments # Nabs = Nexc*(Distance from top of lanex)/(Phosphor thickness)
#where Distance from top of lanex = (SegmentNumber+0.5)/(Nsegments)*(Phosphor thickness)
photonNumberArray.append(photonNumber)
QuantumEfficiency = 0.4
PhotonSignal = np.array(np.multiply(photonNumberArray,QuantumEfficiency))
print(PhotonSignal)
The amount of energy deposited into the CCD depletion region from electrons and x-rays can be read for both the "with vinyl" and "without vinyl" cases:
In [42]:
#With vinyl
ExAbsorbedCCDVinyl = [] #x-rays
ExAbsorbedCCDErrorVinyl = []
EelAbsorbedCCDVinyl = [] #electrons
EelAbsorbedCCDErrorVinyl = []
CCDmass = np.multiply(CCDthickness,1.39740*10**-4) #mass of photoactive layer of CCD
for i in range(len(KEsim)):
ExAbsorbedCCDVinyl.append([])
ExAbsorbedCCDErrorVinyl.append([])
EelAbsorbedCCDVinyl.append([])
EelAbsorbedCCDErrorVinyl.append([])
for j in range(len(CCDthickness)):
#Directory of the input files where vinyl was used in the simulation
if i%2==0: #Energies where vinyl was experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output1/Out_{KE}MeV_{Theta}Degrees_{CCD}umCCD'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)),\
CCD=str(CCDthickness[j]))
else: #Energies where vinyl wasn't experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output2/Out_{KE}MeV_{Theta}Degrees_{CCD}umCCD_inverted'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)),\
CCD=str(CCDthickness[j]))
printflag = False
with open(directory) as searchfile:
firstoccurence = False
for line in searchfile:
left,sep,right = line.partition('cell 7')
if printflag:
Eabs_per_g = float(line[17:28])
Eabs_per_g_Error = float(line[29:35])
if firstoccurence:
ExAbsorbedCCDVinyl[i].append(Eabs_per_g*CCDmass[j])
ExAbsorbedCCDErrorVinyl[i].append(Eabs_per_g_Error*CCDmass[j])
else:
EelAbsorbedCCDVinyl[i].append(Eabs_per_g*CCDmass[j])
EelAbsorbedCCDErrorVinyl[i].append(Eabs_per_g_Error*CCDmass[j])
printflag = False
if sep: # True iff 'cell 6' in line
printflag = True
firstoccurence = not firstoccurence
#Calculate the overall contribution to the signal
Energy_eh_pair = 3.65 #Energy needed to generate electron-hole pair in Si in eV
ElectronSignalVinyl = []
XraySignalVinyl = []
for i in range(len(KEsim)):
ElectronSignalVinyl.append([])
XraySignalVinyl.append([])
ElectronSignalVinyl[i].append(np.multiply(EelAbsorbedCCDVinyl[i],10**6/Energy_eh_pair))
XraySignalVinyl[i].append(np.multiply(ExAbsorbedCCDVinyl[i],10**6/Energy_eh_pair))
ElectronSignalVinyl = np.array(ElectronSignalVinyl)
XraySignalVinyl = np.array(XraySignalVinyl)
TotalSignalVinyl = np.add(ElectronSignalVinyl,XraySignalVinyl)
In [43]:
#Without vinyl
ExAbsorbedCCDNoVinyl = [] #x-rays
ExAbsorbedCCDErrorNoVinyl = []
EelAbsorbedCCDNoVinyl = [] #electrons
EelAbsorbedCCDErrorNoVinyl = []
CCDmass = np.multiply(CCDthickness,1.39740*10**-4) #mass of photoactive layer of CCD
for i in range(len(KEsim)):
ExAbsorbedCCDNoVinyl.append([])
ExAbsorbedCCDErrorNoVinyl.append([])
EelAbsorbedCCDNoVinyl.append([])
EelAbsorbedCCDErrorNoVinyl.append([])
for j in range(len(CCDthickness)):
#Directory of the input files where vinyl was not used in the simulation
if i%2==0: #Energies where vinyl was experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output2/Out_{KE}MeV_{Theta}Degrees_{CCD}umCCD_inverted'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)),\
CCD=str(CCDthickness[j]))
else: #Energies where vinyl wasn't experimentally used
directory = MCNP_Directory + '/MCNP_Decks_Varying_CCD_Thickness/Old_Output/Output1/Out_{KE}MeV_{Theta}Degrees_{CCD}umCCD'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)),\
CCD=str(CCDthickness[j]))
printflag = False
with open(directory) as searchfile:
firstoccurence = False
for line in searchfile:
left,sep,right = line.partition('cell 6')
if printflag:
Eabs_per_g = float(line[17:28])
Eabs_per_g_Error = float(line[29:35])
if firstoccurence:
ExAbsorbedCCDNoVinyl[i].append(Eabs_per_g*CCDmass[j])
ExAbsorbedCCDErrorNoVinyl[i].append(Eabs_per_g_Error*CCDmass[j])
else:
EelAbsorbedCCDNoVinyl[i].append(Eabs_per_g*CCDmass[j])
EelAbsorbedCCDErrorNoVinyl[i].append(Eabs_per_g_Error*CCDmass[j])
printflag = False
if sep: # True iff 'cell 6' in line
printflag = True
firstoccurence = not firstoccurence
#Calculate the overall contribution to the signal
Energy_eh_pair = 3.65 #Energy needed to generate electron-hole pair in Si in eV
ElectronSignalNoVinyl = []
XraySignalNoVinyl = []
for i in range(len(KEsim)):
ElectronSignalNoVinyl.append([])
XraySignalNoVinyl.append([])
ElectronSignalNoVinyl[i].append(np.multiply(EelAbsorbedCCDNoVinyl[i],10**6/Energy_eh_pair))
XraySignalNoVinyl[i].append(np.multiply(ExAbsorbedCCDNoVinyl[i],10**6/Energy_eh_pair))
TotalSignalNoVinyl = np.add(ElectronSignalNoVinyl,XraySignalNoVinyl)
for i in range(len(KEsim)):
TotalSignalNoVinyl[i] += PhotonSignal[i]
Finally, determine the signal ratios
In [44]:
Ratio = [TotalSignalNoVinyl[0][0][0:19]/TotalSignalVinyl[0][0][0:19],\
TotalSignalNoVinyl[1][0][0:19]/TotalSignalVinyl[1][0][0:19],\
TotalSignalNoVinyl[2][0][0:19]/TotalSignalVinyl[2][0][0:19],\
TotalSignalNoVinyl[3][0][0:19]/TotalSignalVinyl[3][0][0:19]]
MCNPSignalRatio = []
for i in range(len(CCDthickness)-1):
MCNPSignalRatio.append([])
for j in range(len(KEsim)):
MCNPSignalRatio[i].append(Ratio[j][i])
and plot the results:
In [46]:
mpl.rcParams.update({'font.size': 24, 'font.family': 'serif'})
sns.set(context='poster',font_scale=1.5)
sns.set_style("darkgrid")
sns.set_palette(palette='deep')
sns.set_color_codes(palette='deep')
plt.figure(figsize=(12,6))
plt.plot(np.multiply(KEsim,1.0),MCNPSignalRatio[0],linewidth=2,label='1 '+ u'\u03bcm')
plt.plot(np.multiply(KEsim,1.0),MCNPSignalRatio[1],linewidth=2,label='2 '+ u'\u03bcm')
plt.plot(np.multiply(KEsim,1.0),MCNPSignalRatio[4],linewidth=2,color='y',label='5 '+ u'\u03bcm')
plt.plot(np.multiply(KEsim,1.0),MCNPSignalRatio[9],linewidth=2,color='c',label='10 '+ u'\u03bcm')
#plt.plot(np.multiply(KEsim,1.0),MCNPSignalRatio[14],linewidth=2,color='purple',label='MCNP 15 '+ u'\u03bcm photoactive region')
plt.plot(np.divide(ChosenEnergies,10**6*q),SignalRatio,linewidth=2,color='r',marker='o',label='Experiment')
plt.xlabel('Electron Energy (keV)')
plt.ylabel('Signal Ratio (No Vinyl / Vinyl)')
plt.xlim(0.33,0.50)
plt.legend(title='Photoactive Region Thickness')
plt.subplots_adjust(left=0.14,bottom=0.15) #Adjust spacing to prevent clipping of x and y labels
#plt.savefig('Signal_Ratio.svg')