In [1]:
from xerawdp_helpers import * # helper functions for retrieving xerawdp data
from Kr83m_Basic import * # pax minitree class for Kr83m data
from cut_helpers import * # functions to apply and plot some event selections
from lce_helpers import * # functions for binning, building map files, and plotting LCE maps
from math import sqrt
from array import array
import numpy as np
import pandas as pd
import glob
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
matplotlib.rc('font', size=16)
plt.rcParams['figure.figsize'] = (20.0,15.0)
import hax
hax.init(main_data_paths=['/home/berget2/scratch-midway/paxProcessed_kr83mDiffusion'])
#hax.ipython.code_hider()
In [ ]:
# from 2015 Kr83m diffusion-mode data
datasets_xerawdp = ['xe100_150413_1839','xe100_150414_1535','xe100_150415_1749',
'xe100_150416_1832','xe100_150417_1648','xe100_150418_1644',
'xe100_150419_0027','xe100_150419_0547','xe100_150419_1611',
'xe100_150420_0304','xe100_150420_1809','xe100_150428_1232',
'xe100_150429_0515','xe100_150429_2203','xe100_150430_1454',
'xe100_150501_0751','xe100_150501_1457','xe100_150502_0802',
'xe100_150503_0110','xe100_150503_1826']
# xerawdp path
xerawdpPath = '/project/lgrandi/tunnell/run_14/NewNN/'
# get xerawdp tree, apply event restrictions, and retrieve desirable data
# see xerawdp_helpers.py for details
xerawdpTree = load_xerawdp_tree(datasets_xerawdp, xerawdpPath)
if 'xerawdp_krRestricted.pkl' not in glob.glob('*'):
df_xerawdp = build_xerawdp_df(xerawdpTree)
df_xerawdp.to_pickle('xerawdp_krRestricted.pkl')
else:
df_xerawdp = pd.read_pickle('xerawdp_krRestricted.pkl')
In [2]:
# datasets processed by pax_4.1.2
datasets_pax = ['xe100_150413_1839','xe100_150414_1535','xe100_150415_1749',
'xe100_150416_1832','xe100_150419_1611','xe100_150420_0304',
'xe100_150420_1809']
for i in range(len(datasets_pax)):
datasets_pax[i]+='_pax4.9.1'
# load minitrees
# throws error when building minitrees for many datasets, I built them one by one
data = hax.minitrees.load(datasets_pax, treemakers=Kr83m_Basic)
df_pax = pd.DataFrame(data[data['s10Time']>=0]) # remove NaNs
In [3]:
# Add some colums to dataframes for easy cut comparison
for df in [df_pax]:
df['s1Dt'] = df['s11Time']-df['s10Time']
df['s1Gap'] = df['s11LeftEdge']-df['s10RightEdge']
df['s20Width'] = df['s20RightEdge']-df['s20LeftEdge']
df['s1sSpan'] = df['s11RightEdge']-df['s10LeftEdge']
df['s1s2Deficit'] = df['s20Width']-df['s1sSpan']
df['s1sRatio'] = df['s11Area']/df['s10Area']
In [4]:
cuts = [ ['s10Coin',2,'none',50,0,100,'channels'],
['s11Coin',2,'none',50,0,100,'channels'],
['s20Area',150,'none',50,0,30000,'PE'],
['s1Gap',0,'none',50,-5000,5000,'ns'],
['s21Area','none',150,50,0,20000,'PE'],
['s1s2Deficit',0,'none',50,-5000,5000,'ns'],
['s1Dt',500,1000,50,-1000,5000,'ns'],
['s1sRatio',0.1,1.0,50,0,1.2,'PE/PE']]
In [5]:
#df_xerawdp_cut = apply_cuts(df_xerawdp,cuts)
#print(len(df_xerawdp_cut.values))
df_pax_cut = apply_cuts(df_pax,cuts)
print(len(df_pax_cut.values))
In [ ]:
#h_dt_xerawdp = ROOT.TH1D('','',70,400,1100)
h_dt_pax = ROOT.TH1D('','',70,400,1100)
#for i in range(len(df_xerawdp_cut.values)):
#h_dt_xerawdp.Fill(df_xerawdp_cut['s1Dt'].values[i])
for i in range(len(df_pax_cut.values)):
h_dt_pax.Fill(df_pax_cut['s1Dt'].values[i])
In [ ]:
def expDec(x,p):
f = p[2]+p[0]*np.exp(-np.log(2)*x[0]/p[1])
return f
In [ ]:
fitFun = ROOT.TF1("fitFun",expDec,500,1000,3)
fitFun.SetParameter(0,200000)
fitFun.SetParameter(1,155)
fitFun.SetParameter(2,10)
c1 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat(0)
c1.SetLogy()
h_dt_pax.GetXaxis().SetTitle('s11Time - s10Time (ns)')
h_dt_pax.GetXaxis().CenterTitle()
h_dt_pax.SetTitle('Pax s1Dt Histogram')
h_dt_pax.SetMaximum(2e4)
h_dt_pax.Draw()
h_dt_pax.Fit("fitFun","","",500,1000)
fit2 = h_dt_pax.GetFunction("fitFun")
chi22 = fit2.GetChisquare()
ndf2 = fit2.GetNDF()
p12 = fit2.GetParameter(1)
e12 = fit2.GetParError(1)
pt2 = ROOT.TPaveText(.58, .68, .88, .88, 'NDC')
pt2.AddText("t_{1/2}=%1.3f"%p12)
pt2.AddText("#sigma=%1.3f"%e12)
pt2.AddText("#chi^{2}/NDF=%1.3f/%i"%(chi22,ndf2))
pt2.Draw()
c1.Print('./KrLce_Figures/dtHists.png')
c1.Clear()
In [ ]:
#plt.rcParams['figure.figsize'] = (20.0, 18.0)
plt.imshow(mpimg.imread('KrLce_Figures/dtHists.png'))
plt.axis('off')
plt.show()
In [ ]:
h_s1_vs_dt = ROOT.TH2D('','',70,500,1000,50,0,120)
for i in range(len(df_pax_cut['s1Dt'].values)):
h_s1_vs_dt.Fill(df_pax_cut['s1Dt'].values[i],df_pax_cut['s11Area'].values[i])
In [ ]:
c2 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat(0)
h_s1_vs_dt.GetXaxis().SetTitle('s1Dt [ns]')
h_s1_vs_dt.GetXaxis().CenterTitle()
h_s1_vs_dt.GetYaxis().SetTitle('s11Area (PE)')
h_s1_vs_dt.GetYaxis().CenterTitle()
h_s1_vs_dt.SetTitle('Pax s11Area vs s1Dt')
h_s1_vs_dt.Draw('colz')
c2.Print('./s1_vs_dt_hist.png')
c2.Clear()
plt.imshow(mpimg.imread('s1_vs_dt_hist.png'))
plt.axis('off')
plt.show()
In [ ]:
for x in range(500, 1000, 25):
hist = h_s1_vs_dt = ROOT.TH1D('','',50,0,100)
for i in range(len(df_pax_cut['s1Dt'].values)):
if df_pax_cut['s1Dt'].values[i] >= x and df_pax_cut['s1Dt'].values[i] < (x+25):
hist.Fill(df_pax_cut['cs11Area'].values[i])
c2 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat('Me')
c2.SetLogy()
hist.GetXaxis().SetTitle('Corrected S1[1] Area [pe]')
hist.GetXaxis().CenterTitle()
hist.SetTitle('Pax cs11Area Histogram: '+str(x)+'ns < dt < '+str(x+25)+'ns')
hist.Draw('hist')
print(x)
c2.Print('./ethans_plots/dthist_'+str(x)+'.png')
c2.Clear()
In [ ]:
from array import array
x = array('d', df_pax_cut['s1Dt'])
y = array('d', df_pax_cut['s11Area'])
g = ROOT.TGraph(len(x), x, y)
c3 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat(0)
c3.SetLogy()
g.GetXaxis().SetTitle('dt [ns]')
g.GetXaxis().CenterTitle()
g.GetYaxis().SetTitle('s11Area [pe]')
g.GetYaxis().CenterTitle()
g.SetTitle('s11Area vs dt')
g.Draw('a*')
c3.Print('./s1_vs_dt.png')
c3.Clear()
plt.imshow(mpimg.imread('s1_vs_dt.png'))
plt.axis('off')
plt.show()
In [ ]:
from math import sqrt
h_s1_vs_r = ROOT.TH2D('','',50,0,25,40,0,120)
count = 0
for i in range(len(df_pax_cut['s1Dt'].values)):
if df_pax_cut['s1Dt'].values[i] <= 550:
count += 1
h_s1_vs_r.Fill( ( sqrt( (df_pax_cut['i0x'].values[i])**2 + (df_pax_cut['i0x'].values[i])**2) ) ,df_pax_cut['s11Area'].values[i])
print (count)
In [ ]:
c5 = ROOT.TCanvas('','',2000,1600)
ROOT.gStyle.SetOptStat(0)
h_s1_vs_r.GetXaxis().SetTitle('r [cm]')
h_s1_vs_r.GetXaxis().CenterTitle()
h_s1_vs_r.GetYaxis().SetTitle('s11Area [pe]')
h_s1_vs_r.GetYaxis().CenterTitle()
h_s1_vs_r.SetTitle('Pax s11Area vs r for dt <= 500 ns')
h_s1_vs_r.Draw('colz')
c5.Print('./s1_vs_r_hist.png')
c5.Clear()
plt.imshow(mpimg.imread('s1_vs_r_hist.png'))
plt.axis('off')
plt.show()
In [ ]:
from math import sqrt
h_s1_vs_r2 = ROOT.TH2D('','',50,0,500,40,0,120)
for i in range(len(df_pax_cut['s1Dt'].values)):
if df_pax_cut['s1Dt'].values[i] <= 550:
h_s1_vs_r2.Fill(((df_pax_cut['i0x'].values[i])**2 + (df_pax_cut['i0x'].values[i])**2),df_pax_cut['s11Area'].values[i])
In [ ]:
c5 = ROOT.TCanvas('','',2000,1600)
ROOT.gStyle.SetOptStat(0)
h_s1_vs_r2.GetXaxis().SetTitle('r^2 [cm^2]')
h_s1_vs_r2.GetXaxis().CenterTitle()
h_s1_vs_r2.GetYaxis().SetTitle('s11Area [pe]')
h_s1_vs_r2.GetYaxis().CenterTitle()
h_s1_vs_r2.GetZaxis().SetRangeUser(0, 500)
h_s1_vs_r2.SetTitle('Pax s11Area vs r^2 for dt <= 500 ns')
h_s1_vs_r2.Draw('colz')
c5.Print('./s1_vs_r2_hist.png')
c5.Clear()
plt.imshow(mpimg.imread('s1_vs_r2_hist.png'))
plt.axis('off')
plt.show()
In [ ]:
h_s1_vs_z = ROOT.TH2D('','',30,-30.3,0,40,0,120)
count = 0
for i in range(len(df_pax_cut['s1Dt'].values)):
if df_pax_cut['s1Dt'].values[i] <= 550:
count +=1
h_s1_vs_z.Fill( df_pax_cut['i0z'].values[i] ,df_pax_cut['cs11Area'].values[i])
print (count)
In [ ]:
c6 = ROOT.TCanvas('','',2000,1600)
ROOT.gStyle.SetOptStat(0)
h_s1_vs_z.GetXaxis().SetTitle('z [cm]')
h_s1_vs_z.GetXaxis().CenterTitle()
h_s1_vs_z.GetYaxis().SetTitle('cs11Area [pe]')
h_s1_vs_z.GetYaxis().CenterTitle()
h_s1_vs_z.GetZaxis().SetRangeUser(0, 500)
h_s1_vs_z.SetTitle('Pax s11Area vs z for dt <= 550 ns')
h_s1_vs_z.Draw('colz')
c6.Print('./s1_vs_z_hist.png')
c6.Clear()
plt.imshow(mpimg.imread('s1_vs_z_hist.png'))
plt.axis('off')
plt.show()
In [6]:
df_newer = df_pax_cut
df_newer = df_newer[ (df_newer['i0z'] > -15) & (df_newer['i0z'] < -1)
& (np.sqrt(df_newer['i0x']**2 + df_newer['i0y']**2) < 13)]
In [35]:
s10hist = ROOT.TH1D('','',100,0,8.5)
for i in range(len(df_newer['cs10Area'].values)):
s10hist.Fill(df_newer['cs10Area'].values[i]/32.1)
In [36]:
c8 = ROOT.TCanvas('','',2000,1600)
ROOT.gStyle.SetOptStat('Me')
s10hist.GetXaxis().SetTitle('cs10 LY [pe/keV]')
s10hist.GetXaxis().CenterTitle()
s10hist.SetTitle('Relative s10 LY Histogram')
s10hist.Draw('hist')
c8.Print('./cs10hist.png')
c8.Clear()
plt.imshow(mpimg.imread('cs10hist.png'))
plt.axis('off')
plt.show()
In [38]:
s11hist = ROOT.TH1D('','',100,0,8.5)
for i in range(len(df_newer['cs11Area'].values)):
s11hist.Fill(df_newer['cs11Area'].values[i]/9.4)
In [39]:
c7 = ROOT.TCanvas('','',2000,1600)
ROOT.gStyle.SetOptStat('Me')
s11hist.GetXaxis().SetTitle('cs11 LY [pe/keV]')
s11hist.GetXaxis().CenterTitle()
s11hist.SetTitle('Relative s11 LY Histogram')
s11hist.Draw('hist')
c7.Print('./cs11hist.png')
c7.Clear()
plt.imshow(mpimg.imread('cs11hist.png'))
plt.axis('off')
plt.show()
In [40]:
f = open('NESTdata.csv')
xvalues = []
yvalues = []
while True:
line = f.readline()
if line == '':
break
line = line.strip()
line = line.split(',')
line[1] = line[1].strip()
xvalues.append(line[0])
yvalues.append(line[1])
f.close()
for i in range(len(xvalues)):
xvalues[i] = float(xvalues[i])
yvalues[i] = float(yvalues[i])
x1 = array('d', xvalues)
y1 = array('d', yvalues)
g1 = ROOT.TGraph(len(x1), x1, y1)
In [41]:
c9 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat(0)
c9.SetLogx()
g1.GetXaxis().SetTitle('beta energy [keV]')
g1.GetXaxis().CenterTitle()
g1.GetYaxis().SetTitle('yield [photons/keV]')
g1.GetYaxis().CenterTitle()
g1.SetTitle('LY as a function of electron energy')
g1.Draw('ac')
c9.Print('./LYenergyNEST.png')
c9.Clear()
plt.imshow(mpimg.imread('LYenergyNEST.png'))
plt.axis('off')
plt.show()
In [14]:
nest_s10_val = g1.Eval(32.1)
our32 = nest_s10_val
print (our32)
scale = nest_s10_val/2.774
print (scale)
In [15]:
nest_s11_val = g1.Eval(9.4)
print (nest_s11_val)
In [16]:
our9 = scale*3.539
print(our9)
In [17]:
newXvalues = [9.4, 32.1]
newYvalues = [our9, our32]
x2 = array('d', newXvalues)
y2 = array('d', newYvalues)
g2 = ROOT.TGraph(len(x2), x2, y2)
g2.SetMarkerColor(2)
In [26]:
c10 = ROOT.TCanvas('','',1600,700)
ROOT.gStyle.SetOptStat(0)
c10.SetLogx()
g1.GetXaxis().SetTitle('beta energy [keV]')
g1.GetXaxis().CenterTitle()
g1.GetYaxis().SetRangeUser(0,54)
g1.GetYaxis().SetTitle('yield [photons/keV]')
g1.GetYaxis().CenterTitle()
g1.SetTitle('LY as a function of electron energy')
g1.Draw('ac')
g2.Draw('*')
c10.Print('./NESTcomp.png')
c10.Clear()
plt.imshow(mpimg.imread('NESTcomp.png'))
plt.axis('off')
plt.show()
In [ ]: