diff --git a/benchmarks/rich/draw_benchmark.py b/benchmarks/rich/draw_benchmark.py
index 4b39934c4c96b1415ff8bebea1ac1a16e06ec237..6f93e4047f2cb9001322c7f3801648177d953fcc 100755
--- a/benchmarks/rich/draw_benchmark.py
+++ b/benchmarks/rich/draw_benchmark.py
@@ -22,7 +22,6 @@ PARTICLE_MASS = {
         'K':  0.49368,
         'p':  0.93827,
         }
-
 # ARGUMENTS
 ################################################################
 
@@ -71,8 +70,8 @@ def make_canv(name, nx, ny):
 canv_dict = {
     "photon_spectra": make_canv("photon_spectra", 2, 2),
     "digitization":   make_canv("digitization",   2, 2),
-    "pidAerogel":     make_canv("pidAerogel",     5, 3),
-    "pidGas":         make_canv("pidGas",         5, 3),
+    "pidAerogel":     make_canv("pidAerogel",     6, 3),
+    "pidGas":         make_canv("pidGas",         6, 3),
     "pidMerged":      make_canv("pidMerged",      2, 3),
 }
 
@@ -93,7 +92,7 @@ def draw_profile(hist, rad, style="BOX"):
     prof = hist.ProfileX("_pfx", 1, -1, "i")
     prof.SetMarkerStyle(r.kFullCircle)
     prof.SetMarkerColor(r.kRed)
-    prof.SetMarkerSize(2)
+    prof.SetMarkerSize(1.5)
     hist_rebinned = hist.Clone(f'{hist.GetName()}_rebin')
     if 'vs_p' in hist.GetName():
         hist_rebinned.RebinX(RADIATORS[rad]['p_rebin'])
@@ -169,48 +168,65 @@ for rad in ["Aerogel", "Gas"]:
     canv.cd(1)
     ana_file.Get(f'{pid_name}/npe_dist_{rad}').Draw()
     canv.cd(2)
-    ana_file.Get(f'{pid_name}/theta_dist_{rad}').Draw()
+    hist = ana_file.Get(f'{pid_name}/theta_dist_{rad}')
+    mean = hist.GetBinCenter(hist.GetMaximumBin())
+    width= hist.GetStdDev()
+    nsigma =3
+    hist.Fit("gaus","","",mean-nsigma*width,mean+nsigma*width)
+    if hist.GetFunction("gaus"):
+        hist.GetFunction("gaus").SetLineWidth(5)        
+        sigma = hist.GetFunction("gaus").GetParameter(2)
+        hist.GetXaxis().SetRangeUser(mean-2*nsigma*sigma,mean+2*nsigma*sigma)
+    hist.Draw()
     canv.cd(3)
-    hist = ana_file.Get(f'{pid_name}/thetaResid_dist_{rad}')
-    hist.DrawCopy()
+    hist = ana_file.Get(f'{pid_name}/speTheta_dist_{rad}')
+    mean = hist.GetBinCenter(hist.GetMaximumBin())
+    width= hist.GetStdDev()
+    print({rad},mean,width)
+    nsigma =3
+    hist.Fit("gaus","","",mean-nsigma*width,mean+nsigma*width)
+    if hist.GetFunction("gaus"):
+        hist.GetFunction("gaus").SetLineWidth(5)                
+        sigma = hist.GetFunction("gaus").GetParameter(2)
+        hist.GetXaxis().SetRangeUser(mean-2*nsigma*sigma,mean+2*nsigma*sigma)
+    hist.Draw()
+    hist.Draw()
     canv.cd(4)
-    hist = ana_file.Get(f'{pid_name}/speResid_dist_{rad}')
+    hist = ana_file.Get(f'{pid_name}/thetaResid_dist_{rad}')
     hist.SetTitle(hist.GetTitle() + " - ZOOM")
     hist.GetXaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
     hist.Draw()
+    nsigma =3
     resid_mode = hist.GetBinCenter(hist.GetMaximumBin())
     resid_dev  = hist.GetStdDev()
-    nsigma = 3
     hist.Fit("gaus", "", "", resid_mode - nsigma*resid_dev, resid_mode + nsigma*resid_dev)
     if hist.GetFunction("gaus"):
         hist.GetFunction("gaus").SetLineWidth(5)
         centralV = hist.GetFunction("gaus").GetParameter(1)
         sigmaV = hist.GetFunction("gaus").GetParameter(2)
-        if centralV > 0.2*sigmaV :
+        if centralV > sigmaV :
             print("dRICH Warning : spe residual plot has starnge central value")
     canv.cd(5)
+    hist = ana_file.Get(f'{pid_name}/speResid_dist_{rad}')
     hist.SetTitle(hist.GetTitle() + " - ZOOM")
     hist.GetXaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
     hist.Draw()
     resid_mode = hist.GetBinCenter(hist.GetMaximumBin())
     resid_dev  = hist.GetStdDev()
-    nsigma = 3
+    nsigma =3
     hist.Fit("gaus", "", "", resid_mode - nsigma*resid_dev, resid_mode + nsigma*resid_dev)
     if hist.GetFunction("gaus"):
         hist.GetFunction("gaus").SetLineWidth(5)
         centralV = hist.GetFunction("gaus").GetParameter(1)
         sigmaV = hist.GetFunction("gaus").GetParameter(2)
-        if centralV > 0.2*sigmaV :
-            print("dRICH Warning : ring residual plot has starnge central value")
+        if centralV > sigmaV :
+            print("dRICH Warning : spe residual plot has starnge central value")
     canv.cd(6)
-    hist = ana_file.Get(f'{pid_name}/speResid_dist_{rad}')
-    hist.DrawCopy()
-    canv.cd(7)
     ana_file.Get(f'{pid_name}/highestWeight_dist_{rad}').Draw()
-    canv.cd(8)
+    canv.cd(7)
     hist = ana_file.Get(f'{pid_name}/npe_vs_p_{rad}')
     draw_profile(hist, rad)
-    canv.cd(9)
+    canv.cd(8)
     hist = ana_file.Get(f'{pid_name}/theta_vs_p_{rad}')
     draw_profile(hist, rad)
     if 'cherenkov_curves' in RADIATORS[rad]:
@@ -224,7 +240,7 @@ for rad in ["Aerogel", "Gas"]:
     hist.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
     draw_profile(hist, rad)
     canv.cd(11)
-    hist = ana_file.Get(f'{pid_name}/thetaResid_vs_p_{rad}')
+    hist = ana_file.Get(f'{pid_name}/speResid_vs_p_{rad}')
     hist.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
     draw_profile(hist, rad)
     canv.cd(12)