Lanex Strip Test

Some quick code for analyzing the data acquired during the Lanex strip test. In this test, the lanex strip used in the electron spectrometer design was partially covered with pieces of electrical tape before being exposed to the electron source. This was done to see how much of the signal is due to Lanex luminesence and how much is due to direct excitation in the CCD by electrons.


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

Espec functions


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)

Data read functions


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)

Load data and subtract background


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

Mesh plot all unsaturated data


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)


/home/drake/.anaconda3/lib/python3.4/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

Plot averaged data


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]:
[<matplotlib.lines.Line2D at 0x7f2873b76e10>]

Interpolate data between regions with and without vinyl


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)

Plot interpolated data


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')


Calculate signal ratios


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)


[  338.   391.   440.   485.   544.   652.   714.   809.   900.  1009.
  1149.  1294.  1428.] 

[ 143.87568911   41.52729299   19.40150658   12.01258033    7.4879801
    5.66289179    6.92871048    4.83261413    6.37340337    5.85732346
    7.96128009    9.35273811   10.56061429]

Determine signal ratios from MCNP output

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]:
<matplotlib.legend.Legend at 0x7f2873d5c1d0>

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)


[ 1560.39655154  2059.36942388  2178.60845554  2124.0356066 ]

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')


/home/drake/.anaconda3/lib/python3.4/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))