diff --git a/code/RDGW150914_ptemcee1.py b/code/RDGW150914_ptemcee1.py
index 21c204f1600b7cd7242529a3db38ee66825659df..ff3e6adc7d40d348c7a5eaf4b16640f5f238dfa8 100755
--- a/code/RDGW150914_ptemcee1.py
+++ b/code/RDGW150914_ptemcee1.py
@@ -2,8 +2,9 @@
 # coding: utf-8
 import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
+
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case varying the fundamental frequencies, tshift = 0. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -41,11 +42,11 @@ tshift=0
 vary_fund = True
 
 #sampler parameters
-npoints=20000 
-nwalkers = 1000
+npoints=200
+nwalkers = 42
 ntemps=16
 ndim = int(4*(nmax+1))
-burnin = 4200  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 42  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
             #This is trivial but often forgotten: this cannot be more than npoints! Usually 1/5~1/4 npoints is what I observe.
 numbins = 42 #corner plot parameter - how many bins you want
 datacolor = '#105670' #'#4fa3a7'
@@ -139,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -161,6 +162,7 @@ def log_likelihood(theta):
         return -np.inf
     return result
 
+
 # Logposterior distribution for the residuals case.
 # The evidence is just a normalization factor
 def log_probability(theta):
@@ -176,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -194,6 +196,7 @@ sampler.run_mcmc(pos,npoints)
 
 
 
+
 #Define labels and start plotting
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
@@ -211,12 +214,12 @@ if vary_fund == False:
 #Chain plot
 fig, axes = plt.subplots(4*(nmax+1), 1, sharex=True, figsize=(12, 9*(nmax+1)))
 for i in range(4*(nmax+1)):
-    axes[i].plot(sampler.chain[0,:, :, i].T, color="k", alpha=0.4)
+    axes[i].plot(sampler.chain[0,:, :, i].T, color="k", alpha=0.4, rasterized=True)
     axes[i].yaxis.set_major_locator(MaxNLocator(5))
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -224,12 +227,20 @@ fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax=
 samples = sampler.chain[0,:, burnin:, :].reshape((-1, ndim))
 #samples for corner plot
 samples_corn = samples if vary_fund == True else np.delete(samples, np.s_[0,2], 1)
+
 #print('Values with peak likelihood:')
 lglk = np.array([log_likelihood(samples[i]) for i in range(len(samples))])
 pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
+if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
+    pk[-2] -= np.pi
+
 #print('pkFalse:')
 #print(pk)
     
@@ -241,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -260,7 +272,7 @@ for yi in range(naxes):
 figcorn.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_corner.pdf', format = 'pdf')
 
 
-#Now export the covariance and correlation matrices, as well as the logevidence
+#Now export the covariance and correlation matrices, as well as the logevidence and autocorrelation time
 #The covariant matrix
 samplesT = samples_corn.T
 cov_matr=np.cov(samplesT)
@@ -279,6 +291,10 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
@@ -318,7 +334,6 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
-
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
 #import os
 #os.system('afplay /System/Library/Sounds/Submarine.aiff')
diff --git a/code/RDGW150914_ptemcee2.py b/code/RDGW150914_ptemcee2.py
index 384fbe97638e504f77dac2d2fdf3b87dc0f09801..546c153955174a94342864fe199c4633a9cf2cd6 100755
--- a/code/RDGW150914_ptemcee2.py
+++ b/code/RDGW150914_ptemcee2.py
@@ -2,8 +2,9 @@
 # coding: utf-8
 import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
+
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case varying the fundamental frequencies, tshift = 19. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -139,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -218,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -226,12 +227,20 @@ fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax=
 samples = sampler.chain[0,:, burnin:, :].reshape((-1, ndim))
 #samples for corner plot
 samples_corn = samples if vary_fund == True else np.delete(samples, np.s_[0,2], 1)
+
 #print('Values with peak likelihood:')
 lglk = np.array([log_likelihood(samples[i]) for i in range(len(samples))])
 pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
+if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
+    pk[-2] -= np.pi
+
 #print('pkFalse:')
 #print(pk)
     
@@ -243,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -281,12 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -320,7 +334,6 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
-
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
 #import os
 #os.system('afplay /System/Library/Sounds/Submarine.aiff')
diff --git a/code/RDGW150914_ptemcee3.py b/code/RDGW150914_ptemcee3.py
index 38b61c8b6df8d2514c96417e00a48c481d3b717a..1e8c2d9acfb5ca713faad1810ea1a229d6398487 100755
--- a/code/RDGW150914_ptemcee3.py
+++ b/code/RDGW150914_ptemcee3.py
@@ -4,7 +4,7 @@ import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
 
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case not varying the fundamental frequencies, tshift = 0. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -42,11 +42,11 @@ tshift=0
 vary_fund = False
 
 #sampler parameters
-npoints=20000 
-nwalkers = 1000
+npoints=200 
+nwalkers = 42
 ntemps=16
 ndim = int(4*(nmax+1))
-burnin = 4200  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 42  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
             #This is trivial but often forgotten: this cannot be more than npoints! Usually 1/5~1/4 npoints is what I observe.
 numbins = 42 #corner plot parameter - how many bins you want
 datacolor = '#105670' #'#4fa3a7'
@@ -140,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -178,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -201,7 +201,7 @@ sampler.run_mcmc(pos,npoints)
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_x = [r'$x_'+str(i)+'$' for i in range (nmax+1)]
-paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_0-\pi$'] + ['$y_'+str(i)+'$' for i in range (1, nmax+1)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_'+str(i)+'-\pi$' for i in range (nmax+1)]
 paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
 #Need to delete alpha_0 and alpha_1 for the corner plot
 paramlabels_corner = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
@@ -219,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -234,9 +234,11 @@ pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
-#y_0 range needs some messaging to make the plot
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
 if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
     samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
     pk[-2] -= np.pi
 
 #print('pkFalse:')
@@ -289,11 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
+
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -327,4 +334,4 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
\ No newline at end of file
diff --git a/code/RDGW150914_ptemcee4.py b/code/RDGW150914_ptemcee4.py
index 8eac5d0295433351bce6303b04ee9b96ae515857..e39042b69fbce490093815ccb504e7b3acdfeb34 100755
--- a/code/RDGW150914_ptemcee4.py
+++ b/code/RDGW150914_ptemcee4.py
@@ -4,7 +4,7 @@ import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
 
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case not varying the fundamental frequencies, tshift = 19. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -140,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -178,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -201,7 +201,7 @@ sampler.run_mcmc(pos,npoints)
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_x = [r'$x_'+str(i)+'$' for i in range (nmax+1)]
-paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_0-\pi$'] + ['$y_'+str(i)+'$' for i in range (1, nmax+1)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_'+str(i)+'-\pi$' for i in range (nmax+1)]
 paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
 #Need to delete alpha_0 and alpha_1 for the corner plot
 paramlabels_corner = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
@@ -219,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -234,9 +234,11 @@ pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
-#y_0 range needs some messaging to make the plot
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
 if vary_fund == False:
-    samples_corn.T[-2] = samples_corn.T[-2] - np.pi
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
     pk[-2] -= np.pi
 
 #print('pkFalse:')
@@ -250,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -288,12 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -327,4 +334,4 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
\ No newline at end of file
diff --git a/code/condor_submit_RdownPtemcee1.sub b/code/condor_submit_RdownPtemcee1.sub
index 52aa79396d33eaae684924e8567d28384073e9b5..d7afc437f7a7bfff552bc471977a01812311e7e4 100755
--- a/code/condor_submit_RdownPtemcee1.sub
+++ b/code/condor_submit_RdownPtemcee1.sub
@@ -13,7 +13,7 @@ notify_user = rl746@cornell.edu
 notification = Complete
 arguments  = "-processid $(Process)" 
 request_memory = 232GB
-request_cpus = 28
+request_cpus = 24
 on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
 accounting_group = aei.dev.test_dynesty
 queue 1
diff --git a/code/condor_submit_RdownPtemcee3.sub b/code/condor_submit_RdownPtemcee3.sub
index 781ec9b69e86d7aa5d8ac400cef03206018bd86a..7e1fbd5644033d1865f25636bac89139a87ffe3b 100755
--- a/code/condor_submit_RdownPtemcee3.sub
+++ b/code/condor_submit_RdownPtemcee3.sub
@@ -13,7 +13,7 @@ notify_user = rl746@cornell.edu
 notification = Complete
 arguments  = "-processid $(Process)" 
 request_memory = 232GB
-request_cpus = 28
+request_cpus = 24
 on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
 accounting_group = aei.dev.test_dynesty
 queue 1