diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 9c6be63475b6e066ece1ceac0aea059b8837a0dc..8f27b36ab8820a66b19133cf8bb526092cd53ad8 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -4,6 +4,7 @@ from GeneratorModules.EvgenAnalysisAlg import EvgenAnalysisAlg
 import ROOT as R
 import numpy as np
 import os
+from math import sqrt
 
 def fix():
     "Python Fixes for HepMC"
@@ -95,22 +96,20 @@ class EvgenValidation(EvgenAnalysisAlg):
             self.hists.add(f"Mass_d{i}", 200, 0, 0.01)            
 
         # Mother
-        self.hists.add("E_M", 100, 0, 10000)
-        self.hists.add("P_M", 100, 0, 10000)
-        self.hists.add("Pz_M", 100, 0, 10000)         
-        self.hists.add("Pt_M", 100, 0, 1)       
-        self.hists.add("Theta_M", 20, 0, 0.001)
+        self.hists.add("E_M", 1000, 0, 10000)
+        self.hists.add("P_M", 1000, 0, 10000)
+        self.hists.add("Pz_M", 1000, 0, 10000)         
+        self.hists.add("Pt_M", 1000, 0, 1)       
+        self.hists.add("Theta_M", 200, 0, 0.02)
         self.hists.add("Phi_M", 16, -3.2, 3.2)
-        self.hists.add("Mass_M", 200, 0, 1)
+        self.hists.add("Mass_M", 200, 0, 2)
         self.hists.add("ThetaVsP_M", arrayX = tbins, arrayY = pbins)
 
         # Vertex
         self.hists.add("Vtx_X", 50, -100, 100)
         self.hists.add("Vtx_Y", 50, -100, 100)
-        # For fluka
-        #self.hists.add("Vtx_X", 100, -3000, 3000)
-        #self.hists.add("Vtx_Y", 100, -3000, 3000)        
-        self.hists.add("Vtx_Z", 50, -1500, 0)
+        self.hists.add("Vtx_Z", 500, -5000, 0)
+        self.hists.add("Vtx_R", 20, 0, 200)        
         self.hists.add("Vtx_XY", 50, -100, 100, 50, -100, 100)        
         
         return StatusCode.Success
@@ -141,7 +140,8 @@ class EvgenValidation(EvgenAnalysisAlg):
         self.hists["Vtx_X"].Fill(v.x(), self.weight)
         self.hists["Vtx_Y"].Fill(v.y(), self.weight)
         self.hists["Vtx_Z"].Fill(v.z(), self.weight)
-        self.hists["Vtx_XY"].Fill(v.x(), v.y(), self.weight)        
+        self.hists["Vtx_XY"].Fill(v.x(), v.y(), self.weight)
+        self.hists["Vtx_R"].Fill(sqrt(v.x()**2 + v.y()**2), self.weight)                
         return
     
 
@@ -184,7 +184,7 @@ if __name__ == "__main__":
     import argparse, sys
     parser = argparse.ArgumentParser(description="Run gen-level validation")
     parser.add_argument("file", nargs="+", help = "full path to imput file")
-    parser.add_argument("--ndaugthers", "-d", default = 2, type = int, help = "Number of daugthers to plot")
+    parser.add_argument("--ndaughters", "-d", default = 2, type = int, help = "Number of daugthers to plot")
     parser.add_argument("--output", "-o",  default = "validation.root", help = "Name of output file")
     parser.add_argument("--mcEventKey", "-k",  default = "BeamTruthEvent", help = "Name of MC collection")
     parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process")    
@@ -218,7 +218,7 @@ if __name__ == "__main__":
     fix()
     
     acc = ComponentAccumulator()
-    valid = EvgenValidation("EvgenValidation", ndaughters =  args.ndaugthers, outname = args.output)
+    valid = EvgenValidation("EvgenValidation", ndaughters =  args.ndaughters, outname = args.output)
     valid.McEventKey = args.mcEventKey
     acc.addEventAlgo(valid)    
     cfg.merge(acc)
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index 8d5a04701c82deedadada27fc76d11fa4a6a8d30..5f702f6953dc891b1878719b36daff4714aee36b 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -14,6 +14,8 @@ def plot(f, name, xtitle, ytitle, xlo = None, xhi = None, ylo = None, yhi = None
 
     if ylo is not None and yhi is not None:
         h.SetAxisRange(ylo, yhi, "Y")
+    elif not logy:
+        h.SetMinimum(0)
 
     if isinstance(r, tuple):
         h.Rebin2D(r[0], r[1])        
@@ -49,7 +51,7 @@ def plot(f, name, xtitle, ytitle, xlo = None, xhi = None, ylo = None, yhi = None
     h.Draw(d)
     return h
 
-def plotn(f, configs, x, y, outname = "valplot"):
+def plotn(f, args, configs, x, y, outname):
 
     c = R.TCanvas()
     c.Divide(x, y)
@@ -62,61 +64,65 @@ def plotn(f, configs, x, y, outname = "valplot"):
         c.cd(i+1)
         c._objs.append(plot(f, *cfg))
         
-    c.Print(f"{outname}.eps")
+    c.Print(f"{args.output}/{outname}.eps")
 
     return
 
 if __name__ == "__main__":
 
+    
     R.gROOT.SetBatch(True)
     R.gStyle.SetOptStat(0)
 
-    fname = "validation.root"
-    f = R.TFile.Open(fname)
-    
-    config = [Hist("P_d0", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
-              Hist("Theta_d0", xtitle = "#theta [rad]", ndiv = -4),
-              Hist("Mass_d0", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4),               
-              Hist("Pt_d0", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
-              Hist("Phi_d0", xtitle = "#phi [rad]"),
-              Hist("ThetaVsP_d0", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz")
-              ]
-
-    plotn(f, config, 3, 2, "daug0")
-
-    config = [Hist("P_d1", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
-              Hist("Theta_d1", xtitle = "#theta [rad]", ndiv = -4),
-              Hist("Mass_d1", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4),               
-              Hist("Pt_d1", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
-              Hist("Phi_d1", xtitle = "#phi [rad]"),
-              Hist("ThetaVsP_d1", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz")
-              ]
-
-    plotn(f, config, 3, 2, "daug1")    
-
-    config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
-              Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4),
-              Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.05),               
-              Hist("Pt_M", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
-              Hist("Phi_M", xtitle = "#phi [rad]"),
-              Hist("ThetaVsP_M", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz")
+    import argparse, sys, os
+    parser = argparse.ArgumentParser(description="Run gen-level validation plotting")
+    parser.add_argument("file", help = "full path to imput file")
+    parser.add_argument("--output", "-o",  default = "valplot", help = "Name of output directory")
+    parser.add_argument("--ndaughters", "-d", default = 2, type = int, help = "Number of daugthers to plot")
+    args = parser.parse_args()
+
+    if not os.path.exists(args.output):
+        os.mkdir(args.output)
+
+    print (args.file)
+    f = R.TFile.Open(args.file)
+
+    for i in range(args.ndaughters):
+        config = [Hist(f"P_d{i}", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
+                  Hist(f"Theta_d{i}", xtitle = "#theta [rad]", ndiv = -4),
+                  Hist(f"Mass_d{i}", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4),               
+                  Hist(f"Pt_d{i}", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
+                  Hist(f"Phi_d{i}", xtitle = "#phi [rad]"),
+                  Hist(f"ThetaVsP_d{i}", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz")
+                  ]
+
+        plotn(f, args, config, 3, 2, f"daug{i}")
+        
+    config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, xlo = 0, xhi = 2000),
+              Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4, r = 10),
+              Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.2, ndiv = 5),               
+              Hist("Pt_M", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 50),
+              Hist("Phi_M", xtitle = "#phi [rad]", r = 2),
+              Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz")
               ]
 
-    plotn(f, config, 3, 2, "mother")
+    plotn(f, args, config, 3, 2, "mother")
 
-    plotn(f, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid") 
+    plotn(f, args, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid") 
 
-    config = [Hist("ThetaVsP_M", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
-              Hist("ThetaVsP_d0", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
-              Hist("ThetaVsP_d1", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz")
-              ]                   
-
-    plotn(f, config, 2, 2, "twod")
+    
+#     config = [Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
+#               Hist("ThetaVsP_d0", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
+#               Hist("ThetaVsP_d1", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz")
+#               ]                   
+# 
+#     plotn(f, args, config, 2, 2, "twod")
 
     config = [Hist("Vtx_X", xtitle = "x [mm]", r = 5),
               Hist("Vtx_Y", xtitle = "y [mm]", r = 5),
-              Hist("Vtx_Z", xtitle = "z [mm]", r = 5),
-              Hist("Vtx_XY", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5))
+              Hist("Vtx_Z", xtitle = "z [mm]", r = 5, ndiv = 5),
+              Hist("Vtx_XY", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
+              Hist("Vtx_R", xtitle = "r [mm]", r = 5, ndiv = 5)
               ]
 
-    plotn(f, config, 2, 2, "vtx")
+    plotn(f, args, config, 3, 2, "vtx")