diff --git a/gate-pbt/analysis/analysis.py b/gate-pbt/analysis/analysis.py
index f916d89..7c9a4a1 100644
--- a/gate-pbt/analysis/analysis.py
+++ b/gate-pbt/analysis/analysis.py
@@ -25,6 +25,7 @@
 import dicomtomhd
 import gamma
 import overrides
+from uncertainty_stats import getUncertaintyPlots
 
 
 
@@ -153,9 +154,10 @@ def full_analysis( outputdir ):
     
     
     print("Fields found: ", fieldnames)
-    
     for field in fieldnames:
-        
+
+      
+            
         print("\nAnalyzing field: ", field)
 
         print("  Merging results...")
@@ -181,7 +183,7 @@ def full_analysis( outputdir ):
 
         
         path_to_dcmdose = mhdtodicom.get_dcm_file_path( outputdir, beamref )
-        tps_dose = dicomtomhd.dcm2mhd( path_to_dcmdose )        
+        tps_dose = dicomtomhd.dcm2mhd( path_to_dcmdose )      
         print("  Overriding TPS outside of zSurface to zero for gamma analysis")
         struct_file = mhdtodicom.get_struct_file_path( outputdir )
         tps_dose = overrides.set_external_dose_zero( tps_dose, struct_file, "zSurface" )    ## hard-coded
@@ -224,7 +226,8 @@ def full_analysis( outputdir ):
             #itk.imwrite(dose_none_ext, path_to_none_ext)      
             #mhdtodicom.mhd2dcm(dose_none_ext, path_to_dcmdose, join(outputdir, field+"_DoseToWater_NoneExt.dcm") )
    
-               
+            #'''
+            
             print("  Performing gamma analysis 3%/3mm: post-sim D2W vs Eclipse")
             gamma_img = gamma.gamma_image(  d2wimg, tps_dose, 3, 3 )
             itk.imwrite(gamma_img, join(outputdir, field+"_Gamma_33.mhd") )
@@ -242,7 +245,8 @@ def full_analysis( outputdir ):
             itk.imwrite(gamma_img_22, join(outputdir, field+"_Gamma_22.mhd") )
             pass_rate = gamma.get_pass_rate( gamma_img_22 )
             print("   *** Gamma pass rate @ 2%/2mm = {}%".format( round(pass_rate,2) ))
-            
+
+            #'''
 
         
         dose2water = field+"_merged-DoseToWater.mhd"
@@ -250,6 +254,7 @@ def full_analysis( outputdir ):
             print("\n")
             print("  Scaling merged-DoseToWater.mhd")
             doseimg_path = join(outputdir, dose2water)
+
             
             # Read DoseMask to set dose outside zSurface to zero
             dosemask = itk.imread(join(parentdir,"data","DoseMask.mhd"))  # TODO; read from config file
@@ -272,6 +277,24 @@ def full_analysis( outputdir ):
             print("   *** Gamma GD2W pass rate @ 3%/3mm = {}%".format( round(pass_rate,2) ))          
             print("  Converting gamma image to dicom")
             gamma_dcm = join(outputdir, field+"_Gamma_GD2W_33.dcm")
+        
+
+        # Maybe we add uncertainy metrics here?
+
+        dose = field+"_merged-Uncertainty.mhd"
+        if dose in [basename(f) for f in mergedfiles]:
+            print("\n")
+            print("Calculating uncertainty")
+            doseUncert_path = join(outputdir, dose)
+            dose2W = d2wimg
+            print("Number of Primaries Simulated", nsim)
+            results = getUncertaintyPlots(doseUncert_path, dose2W, nsim)
+            
+
+        # Dose To Water MHD file
+        
+
+
 
 
         
diff --git a/gate-pbt/analysis/gamma.py b/gate-pbt/analysis/gamma.py
index 89f30cf..ba0a86b 100644
--- a/gate-pbt/analysis/gamma.py
+++ b/gate-pbt/analysis/gamma.py
@@ -20,7 +20,7 @@
 import itk
 
 #from gatetools import gamma_index as gi
-import gamma_index
+import gamma_index as gamma_index
 
 import reorientate
 
@@ -39,7 +39,7 @@ def gamma_image( ref_dose, target_dose, dta_val, dd_val ):
     Accepts ITK-like image, or path to image
     Returns image matching dimensions of target
     Set ref -> MC dose, target -> TPS dose
-    """    
+    """   
     ref, targ = None, None 
     if type(ref_dose)==str:
         #Assume we have file path
diff --git a/gate-pbt/analysis/gamma_index.py b/gate-pbt/analysis/gamma_index.py
index 6c79327..8fd5925 100644
--- a/gate-pbt/analysis/gamma_index.py
+++ b/gate-pbt/analysis/gamma_index.py
@@ -17,7 +17,6 @@
 import numpy as np
 import itk
 from scipy.spatial import cKDTree
-from scipy.interpolate import RegularGridInterpolator
 
 def get_gamma_index(ref,target,**kwargs):
     """
@@ -99,7 +98,7 @@ def GetGamma(d0, d1, x0, x1, y0, y1, z0, z1, Max, dd, dta, gamma_method):
     else:
         print('Please specify gamma method correctly, local or global.')
         return False
-
+    
     return np.sqrt(
         (d1 - d0) ** 2 / (0.01 * dd * norm_val) ** 2 +
         (x1 - x0) ** 2 / dta ** 2 +
@@ -181,7 +180,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
     
 
     
-    gamma, total = 0,0 # for gamma calc in testing
+    gamma, gammaLocal, total = 0, 0, 0 # for gamma calc in testing
     
     # Loop over each target position 
     for Xt in x_tgt:
@@ -214,11 +213,13 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
 
                
                 found = False  # Flag to indicate when to exit all loops
+                foundLoc = False
 
                 # iterate over 4 closest references voxels
                 for idx in indices:
-                    if found:
-                        break  # Exit the outer loop if found
+                    if foundLoc:
+                        if found:
+                            break  # Exit the outer loop if found
                     
                     closest_ref_point = ref_coords[idx]  # (x, y, z) coordinates of closest reference point
                     X, Y, Z = closest_ref_point[0], closest_ref_point[1], closest_ref_point[2]
@@ -235,13 +236,21 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
                     # strictly speaking, this is not a first filter, as we only loop over 4 voxels
                     # in reality, we would need to loop over all voxels within dta and check for filter 1
                     Gamma = GetGamma(tgt_value, ref_value, Xt, X, Yt, Y, Zt, Z, max_targ, dd, dta, gamma_method)
+                    GammaLocal = GetGamma(tgt_value, ref_value, Xt, X, Yt, Y, Zt, Z, max_targ, dd, dta, 'local')
+                    if GammaLocal <=1:
+                        gammaLocal +=1
+                        foundLoc = True
+                        
                     if Gamma <= 1: 
                         gamma += 1
-                        found = True # we have passed this target voxel
+                        found = True
                         # update gamma array values for each method
                         gamma_array[tuple(tgt_index)] = Gamma
-                        break  # Exit the outer loop
+                        #break  # Exit the outer loop
                     
+                    if foundLoc:
+                        if found:
+                            break
 
                     dirns = [+1, -1]
 
@@ -251,6 +260,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
                         # loop over both directions, i.e. +X or -X
                         for dirn in dirns:
                             x,y,z = X,Y,Z
+                            xLoc,yLoc,zLoc = X,Y,Z
                             # interpolating in X
                             if D == 0:
                                 # check we are not at a boundary
@@ -263,6 +273,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
 
                                 #calculate the optimum x position, and the corresponding dose value
                                 x, interpolated_value = min_gamma(X, X + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Xt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method)
+                                xLoc, interpolated_valueLoc = min_gamma(X, X + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Xt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local')
 
                                 # simple check to see if x is within bounds of approximation
                                 if x == False:
@@ -275,6 +286,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
                                     continue
                                 ref_value_next = referenceArray[i_ref, int(j_ref + dirn), k_ref]
                                 y, interpolated_value = min_gamma(Y, Y + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Yt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method)
+                                yLoc, interpolated_valueLoc = min_gamma(Y, Y + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Yt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local')
                                 if y == False:
                                     continue
                             # interpolating in Z
@@ -287,25 +299,36 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0
                                 if ref_value == ref_value_next:
                                     continue
                                 z, interpolated_value = min_gamma(Z, Z + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Zt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method)
+                                zLoc, interpolated_valueLoc = min_gamma(Z, Z + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Zt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local')
                                 if z == False:
                                     continue
 
                             # simple gamma calc function
                             Gamma = GetGamma(tgt_value, interpolated_value, Xt, x, Yt, y, Zt, z, max_targ, dd, dta, gamma_method)
-
-                            # do they pass?
+                            GammaLocal = GetGamma(tgt_value, interpolated_valueLoc, Xt, xLoc, Yt, yLoc, Zt, zLoc, max_targ, dd, dta, 'local')
+                            
                             
+                            
+                            # do they pass?
+                            if GammaLocal <=1:
+                                gammaLocal +=1
+                                foundLoc =True
                             if Gamma <= 1:
                                 gamma += 1
                                 found = True
                                 gamma_array[tuple(tgt_index)] = Gamma
-                                break  # Exit the inner loop if condition is met
+                            if foundLoc:
+                                if found:
+                                    break
+                        if foundLoc:
+                            if found:
+                                break  # Exit the outer loop if inner loop condition was met
+                    if foundLoc:
                         if found:
-                            break  # Exit the outer loop if inner loop condition was met
-                    if found:
-                        break
+                            break
     
-    print('Linear optimum  linear interpolation method = ',100.0 - 100*gamma_array[gamma_array>1].size / gamma_array[gamma_array>0].size, '%')
+    print('Global Gamma - linear interpolation method = ',100.0 - 100*gamma_array[gamma_array>1].size / gamma_array[gamma_array>0].size, '%')
+    print('Local Gamma - linear interpolation method = ', gammaLocal/total * 100, '%')
     
     # convert to image, note we copy information from targetimage
     gimg=itk.image_from_array(gamma_array.swapaxes(0,2).astype(np.float32).copy())
@@ -356,7 +379,7 @@ def min_gamma(x1, x2, y1, y2, tx, ty, dd, dta, Max, spacing, gamma_method):
     # optimal x position for local or global gamma
    
     if gamma_method == 'local':
-        x_opto = (grad*ty/(dd*0.01*ty)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*ty)**2)/(grad**2/(dd*0.01*ty)**2 + 1/dta**2)
+        x_opto= (grad*ty/(dd*0.01*ty)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*ty)**2)/(grad**2/(dd*0.01*ty)**2 + 1/dta**2)
     elif gamma_method == 'global':
         x_opto = (grad*ty/(dd*0.01*Max)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*Max)**2)/(grad**2/(dd*0.01*Max)**2 + 1/dta**2)
     else:
diff --git a/gate-pbt/analysis/mergeresults.py b/gate-pbt/analysis/mergeresults.py
index c173d67..7293c5f 100644
--- a/gate-pbt/analysis/mergeresults.py
+++ b/gate-pbt/analysis/mergeresults.py
@@ -84,7 +84,7 @@ def combine_uncertainty(dosefiles, dosesquaredfiles, statfiles, output):
         #treatment planning" and dividing by dose/N for relative uncertainty...
         if sumdosesq[i]!=0 and sumdose[i]!=0 and N>1: 
             uncertainty[i] = math.sqrt( 1.0/(N-1) * (sumdosesq[i]/N - (sumdose[i]/N)**2) ) / (sumdose[i]/N)
-
+    
     combinedUncert = uncertainty.reshape( shape )
     uncertimg = itk.image_from_array( combinedUncert.astype(np.float32)  ) ## ITK CANNOT WRITE DOUBLES, MUST CAST TO FLOAT
     uncertimg.CopyInformation( img1 )
diff --git a/gate-pbt/analysis/mhdtodicom.py b/gate-pbt/analysis/mhdtodicom.py
index 050a1a8..9d57ab9 100644
--- a/gate-pbt/analysis/mhdtodicom.py
+++ b/gate-pbt/analysis/mhdtodicom.py
@@ -109,6 +109,7 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None):
         print("  Dose scaling of {} used in mhd2dcm".format(dosescaling))
     
     dcm = pydicom.dcmread(dcmFile)
+    
     mhd=None
     if type(mhdFile)==str:
         # Assume file path
@@ -159,7 +160,7 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None):
     # TODO: check this is safe
     # TODO: should this ever negative? - NO!
     dcm.GridFrameOffsetVector = [ x*mhd.GetSpacing()[2] for x in range(mhdpix.shape[0]) ]
-         
+
     dose_abs = mhdpix * dosescaling
     
     scale_to_int = 1E4   # Dicom stores array of integers only
@@ -174,5 +175,4 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None):
     #dosescaling = dcm.DoseGridScaling
     #FrameIncrementPointer ?
    
-    dcm.save_as( output )
-
+    dcm.save_as( output)
diff --git a/gate-pbt/analysis/resampled_dosemask.mhd b/gate-pbt/analysis/resampled_dosemask.mhd
new file mode 100644
index 0000000..83a40d5
--- /dev/null
+++ b/gate-pbt/analysis/resampled_dosemask.mhd
@@ -0,0 +1,13 @@
+ObjectType = Image
+NDims = 3
+BinaryData = True
+BinaryDataByteOrderMSB = False
+CompressedData = False
+TransformMatrix = 1 0 0 0 1 0 0 0 1
+Offset = -102.29000000000001 33.5655 -44
+CenterOfRotation = 0 0 0
+AnatomicalOrientation = RAI
+ElementSpacing = 1.96234 1.97088 2
+DimSize = 106 55 85
+ElementType = MET_UCHAR
+ElementDataFile = resampled_dosemask.raw
diff --git a/gate-pbt/analysis/resampled_dosemask.raw b/gate-pbt/analysis/resampled_dosemask.raw
new file mode 100644
index 0000000..9e6d9cd
Binary files /dev/null and b/gate-pbt/analysis/resampled_dosemask.raw differ
diff --git a/gate-pbt/analysis/uncertainty_stats.py b/gate-pbt/analysis/uncertainty_stats.py
index c08b8a4..afa17aa 100644
--- a/gate-pbt/analysis/uncertainty_stats.py
+++ b/gate-pbt/analysis/uncertainty_stats.py
@@ -9,7 +9,7 @@
 import numpy as np
 import pydicom
 import matplotlib.pyplot as plt
-
+import os
 import gatetools.roi_utils as roiutils
 
 
@@ -61,10 +61,184 @@ def main():
     get_stats_struct(uncertimgpath, rsdcmfile, STRUCT_NAME)
 
 
+import numpy as np
+import matplotlib.pyplot as plt
+import itk
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import itk
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import itk
+
+import numpy as np
+import matplotlib.pyplot as plt
+import itk
+
+def getUncertaintyPlots(uncertainty_path, dose_path, nsim):
+    """
+    Analyze and visualize dose uncertainty data, focusing on clinically relevant regions.
+    """
+    try:
+        ImageType = itk.Image[itk.F, 3]
+        
+        # Read uncertainty image
+        uncert_reader = itk.ImageFileReader[ImageType].New()
+        uncert_reader.SetFileName(uncertainty_path)
+        uncert_reader.Update()
+        doseUncert_img = uncert_reader.GetOutput()
+        
+        # Read dose image
+        dose_reader = itk.ImageFileReader[ImageType].New()
+        dose_reader.SetFileName(dose_path)
+        dose_reader.Update()
+        dose_img = dose_reader.GetOutput()
+        
+    except RuntimeError as e:
+        print(f"Error loading images: {e}")
+        raise
+    
+    # Convert to numpy arrays
+    doseUncert = itk.array_view_from_image(doseUncert_img).swapaxes(0, 2)
+    dose = itk.array_view_from_image(dose_img).swapaxes(0, 2)
+
+
+
+    UncertOrigin = np.array(doseUncert_img.GetOrigin())  # Adjusted for swapped axes
+    UncertSpacing = np.array(doseUncert_img.GetSpacing())
+    doseOrigin = np.array(dose_img.GetOrigin())
+    doseSpacing = np.array(dose_img.GetSpacing())
+
+    '''
+    
+    # Find maximum dose for threshold calculations
+    max_dose = np.max(dose)
+    
+    # Create masks for different dose levels
+    mask_1pct = dose > (0.01 * max_dose)
+    mask_10pct = dose > (0.10 * max_dose)
+    mask_50pct = dose > (0.50 * max_dose)
+    
+    # Calculate statistics for different regions
+    stats = {
+        'all_voxels': {
+            'mean': np.mean(doseUncert),
+            'median': np.median(doseUncert),
+            'std': np.std(doseUncert)
+        },
+        'above_1pct': {
+            'mean': np.mean(doseUncert[mask_1pct]),
+            'median': np.median(doseUncert[mask_1pct]),
+            'std': np.std(doseUncert[mask_1pct])
+        },
+        'above_10pct': {
+            'mean': np.mean(doseUncert[mask_10pct]),
+            'median': np.median(doseUncert[mask_10pct]),
+            'std': np.std(doseUncert[mask_10pct])
+        },
+        'above_50pct': {
+            'mean': np.mean(doseUncert[mask_50pct]),
+            'median': np.median(doseUncert[mask_50pct]),
+            'std': np.std(doseUncert[mask_50pct])
+        }
+    }
+    
+    # Create figure with multiple plots
+    fig = plt.figure(figsize=(15, 15))
+    
+    # Adjust the layout to prevent overlap
+    plt.subplots_adjust(hspace=0.4, wspace=0.3)
+    
+    # Plot 1: All voxels
+    ax1 = plt.subplot(221)
+    ax1.hist(doseUncert.flatten(), bins=100, density=True, alpha=0.75, color='skyblue')
+    ax1.set_title('All Voxels', pad=20, fontsize=12, fontweight='bold')
+    ax1.set_xlabel('Relative Uncertainty', fontsize=10)
+    ax1.set_ylabel('Density', fontsize=10)
+    mean_line = ax1.axvline(stats['all_voxels']['mean'], color='red', linestyle='--', 
+                label=f"Mean: {stats['all_voxels']['mean']:.3f}")
+    ax1.legend(loc='upper right')
+    
+    # Plot 2: >1% max dose
+    ax2 = plt.subplot(222)
+    ax2.hist(doseUncert[mask_1pct].flatten(), bins=100, density=True, alpha=0.75, color='lightgreen')
+    ax2.set_title('Voxels >1% of Max Dose', pad=20, fontsize=12, fontweight='bold')
+    ax2.set_xlabel('Relative Uncertainty', fontsize=10)
+    ax2.set_ylabel('Density', fontsize=10)
+    mean_line = ax2.axvline(stats['above_1pct']['mean'], color='red', linestyle='--', 
+                label=f"Mean: {stats['above_1pct']['mean']:.3f}")
+    ax2.legend(loc='upper right')
+    
+    # Plot 3: >10% max dose
+    ax3 = plt.subplot(223)
+    ax3.hist(doseUncert[mask_10pct].flatten(), bins=100, density=True, alpha=0.75, color='salmon')
+    ax3.set_title('Voxels >10% of Max Dose', pad=20, fontsize=12, fontweight='bold')
+    ax3.set_xlabel('Relative Uncertainty', fontsize=10)
+    ax3.set_ylabel('Density', fontsize=10)
+    mean_line = ax3.axvline(stats['above_10pct']['mean'], color='red', linestyle='--', 
+                label=f"Mean: {stats['above_10pct']['mean']:.3f}")
+    ax3.legend(loc='upper right')
+    
+    # Plot 4: >50% max dose
+    ax4 = plt.subplot(224)
+    ax4.hist(doseUncert[mask_50pct].flatten(), bins=100, density=True, alpha=0.75, color='plum')
+    ax4.set_title('Voxels >50% of Max Dose', pad=20, fontsize=12, fontweight='bold')
+    ax4.set_xlabel('Relative Uncertainty', fontsize=10)
+    ax4.set_ylabel('Density', fontsize=10)
+    mean_line = ax4.axvline(stats['above_50pct']['mean'], color='red', linestyle='--', 
+                label=f"Mean: {stats['above_50pct']['mean']:.3f}")
+    ax4.legend(loc='upper right')
+    
+    # Add grid to all plots
+    for ax in [ax1, ax2, ax3, ax4]:
+        ax.grid(True, alpha=0.3)
+    
+    # Add a main title to the figure
+    plt.suptitle('Dose Uncertainty Distribution Analysis', fontsize=14, fontweight='bold', y=0.95)
+    
+    plt.show()
+    '''
+
+    top_indices = np.argsort(dose.flatten())[-10:]
+    max_dose_median = np.median(dose.flatten()[top_indices])
+
+    fifty_mask = dose >= 0.5 * (max_dose_median)
+    fifty_uncert = doseUncert[fifty_mask] * 100
+    
+    ninety_mask = dose >= 0.9 * (max_dose_median)
+    ninety_uncert = doseUncert[ninety_mask] * 100
+    
+
+    
+    results = {
+        'top50%_mean': np.mean(fifty_uncert),
+        'top50%_median': np.median(fifty_uncert),
+        'top90%_mean': np.mean(ninety_uncert),
+        'top90%_median': np.median(ninety_uncert),
+        'voxel_size_x': doseSpacing[0],
+        'voxel_size_y': doseSpacing[1],
+        'voxel_size_z': doseSpacing[2],
+        'n_primaries': nsim,
+        'total_voxels': dose.size
+    }
+    
+    output_path = os.path.join(os.getcwd(), 'uncertainty_metrics.txt')
+    with open(output_path, 'w') as f:
+        headers = '\t'.join(results.keys())
+        values = '\t'.join(str(x) for x in results.values())
+        f.write(headers + '\n')
+        f.write(values + '\n')
+    
+    print(f"File saved to: {output_path}")
+    return results
 
 
-if __name__=="__main__":
-    main()
+#if __name__=="__main__":
+#    main()