Skip to content
Snippets Groups Projects
Commit 51c5ec1d authored by Carl Gwilliam's avatar Carl Gwilliam
Browse files

Update validation scripts

parent 85979230
No related branches found
No related tags found
3 merge requests!244Signal generation,!243Signal generation,!242MDC Updates
Hide whitespace changes
Inline Side-by-side
Generators/ForeseeGenerator/python/Validate.py
+ 13
− 13
View file @ 51c5ec1d
...@@ -4,6 +4,7 @@ from GeneratorModules.EvgenAnalysisAlg import EvgenAnalysisAlg ...@@ -4,6 +4,7 @@ from GeneratorModules.EvgenAnalysisAlg import EvgenAnalysisAlg
import ROOT as R import ROOT as R
import numpy as np import numpy as np
import os import os
from math import sqrt
def fix(): def fix():
"Python Fixes for HepMC" "Python Fixes for HepMC"
...@@ -95,22 +96,20 @@ class EvgenValidation(EvgenAnalysisAlg): ...@@ -95,22 +96,20 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists.add(f"Mass_d{i}", 200, 0, 0.01) self.hists.add(f"Mass_d{i}", 200, 0, 0.01)
# Mother # Mother
self.hists.add("E_M", 100, 0, 10000) self.hists.add("E_M", 1000, 0, 10000)
self.hists.add("P_M", 100, 0, 10000) self.hists.add("P_M", 1000, 0, 10000)
self.hists.add("Pz_M", 100, 0, 10000) self.hists.add("Pz_M", 1000, 0, 10000)
self.hists.add("Pt_M", 100, 0, 1) self.hists.add("Pt_M", 1000, 0, 1)
self.hists.add("Theta_M", 20, 0, 0.001) self.hists.add("Theta_M", 200, 0, 0.02)
self.hists.add("Phi_M", 16, -3.2, 3.2) 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) self.hists.add("ThetaVsP_M", arrayX = tbins, arrayY = pbins)
# Vertex # Vertex
self.hists.add("Vtx_X", 50, -100, 100) self.hists.add("Vtx_X", 50, -100, 100)
self.hists.add("Vtx_Y", 50, -100, 100) self.hists.add("Vtx_Y", 50, -100, 100)
# For fluka self.hists.add("Vtx_Z", 500, -5000, 0)
#self.hists.add("Vtx_X", 100, -3000, 3000) self.hists.add("Vtx_R", 20, 0, 200)
#self.hists.add("Vtx_Y", 100, -3000, 3000)
self.hists.add("Vtx_Z", 50, -1500, 0)
self.hists.add("Vtx_XY", 50, -100, 100, 50, -100, 100) self.hists.add("Vtx_XY", 50, -100, 100, 50, -100, 100)
return StatusCode.Success return StatusCode.Success
...@@ -141,7 +140,8 @@ class EvgenValidation(EvgenAnalysisAlg): ...@@ -141,7 +140,8 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists["Vtx_X"].Fill(v.x(), self.weight) self.hists["Vtx_X"].Fill(v.x(), self.weight)
self.hists["Vtx_Y"].Fill(v.y(), self.weight) self.hists["Vtx_Y"].Fill(v.y(), self.weight)
self.hists["Vtx_Z"].Fill(v.z(), 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 return
...@@ -184,7 +184,7 @@ if __name__ == "__main__": ...@@ -184,7 +184,7 @@ if __name__ == "__main__":
import argparse, sys import argparse, sys
parser = argparse.ArgumentParser(description="Run gen-level validation") parser = argparse.ArgumentParser(description="Run gen-level validation")
parser.add_argument("file", nargs="+", help = "full path to imput file") 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("--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("--mcEventKey", "-k", default = "BeamTruthEvent", help = "Name of MC collection")
parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process") parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process")
...@@ -218,7 +218,7 @@ if __name__ == "__main__": ...@@ -218,7 +218,7 @@ if __name__ == "__main__":
fix() fix()
acc = ComponentAccumulator() acc = ComponentAccumulator()
valid = EvgenValidation("EvgenValidation", ndaughters = args.ndaugthers, outname = args.output) valid = EvgenValidation("EvgenValidation", ndaughters = args.ndaughters, outname = args.output)
valid.McEventKey = args.mcEventKey valid.McEventKey = args.mcEventKey
acc.addEventAlgo(valid) acc.addEventAlgo(valid)
cfg.merge(acc) cfg.merge(acc)
......
Generators/ForeseeGenerator/share/plot_validation.py
+ 48
− 42
View file @ 51c5ec1d
...@@ -14,6 +14,8 @@ def plot(f, name, xtitle, ytitle, xlo = None, xhi = None, ylo = None, yhi = None ...@@ -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: if ylo is not None and yhi is not None:
h.SetAxisRange(ylo, yhi, "Y") h.SetAxisRange(ylo, yhi, "Y")
elif not logy:
h.SetMinimum(0)
if isinstance(r, tuple): if isinstance(r, tuple):
h.Rebin2D(r[0], r[1]) h.Rebin2D(r[0], r[1])
...@@ -49,7 +51,7 @@ def plot(f, name, xtitle, ytitle, xlo = None, xhi = None, ylo = None, yhi = None ...@@ -49,7 +51,7 @@ def plot(f, name, xtitle, ytitle, xlo = None, xhi = None, ylo = None, yhi = None
h.Draw(d) h.Draw(d)
return h return h
def plotn(f, configs, x, y, outname = "valplot"): def plotn(f, args, configs, x, y, outname):
c = R.TCanvas() c = R.TCanvas()
c.Divide(x, y) c.Divide(x, y)
...@@ -62,61 +64,65 @@ def plotn(f, configs, x, y, outname = "valplot"): ...@@ -62,61 +64,65 @@ def plotn(f, configs, x, y, outname = "valplot"):
c.cd(i+1) c.cd(i+1)
c._objs.append(plot(f, *cfg)) c._objs.append(plot(f, *cfg))
c.Print(f"{outname}.eps") c.Print(f"{args.output}/{outname}.eps")
return return
if __name__ == "__main__": if __name__ == "__main__":
R.gROOT.SetBatch(True) R.gROOT.SetBatch(True)
R.gStyle.SetOptStat(0) R.gStyle.SetOptStat(0)
fname = "validation.root" import argparse, sys, os
f = R.TFile.Open(fname) parser = argparse.ArgumentParser(description="Run gen-level validation plotting")
parser.add_argument("file", help = "full path to imput file")
config = [Hist("P_d0", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5), parser.add_argument("--output", "-o", default = "valplot", help = "Name of output directory")
Hist("Theta_d0", xtitle = "#theta [rad]", ndiv = -4), parser.add_argument("--ndaughters", "-d", default = 2, type = int, help = "Number of daugthers to plot")
Hist("Mass_d0", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4), args = parser.parse_args()
Hist("Pt_d0", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
Hist("Phi_d0", xtitle = "#phi [rad]"), if not os.path.exists(args.output):
Hist("ThetaVsP_d0", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz") os.mkdir(args.output)
]
print (args.file)
plotn(f, config, 3, 2, "daug0") f = R.TFile.Open(args.file)
config = [Hist("P_d1", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5), for i in range(args.ndaughters):
Hist("Theta_d1", xtitle = "#theta [rad]", ndiv = -4), config = [Hist(f"P_d{i}", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
Hist("Mass_d1", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4), Hist(f"Theta_d{i}", xtitle = "#theta [rad]", ndiv = -4),
Hist("Pt_d1", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5), Hist(f"Mass_d{i}", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.001, ndiv = 4),
Hist("Phi_d1", xtitle = "#phi [rad]"), Hist(f"Pt_d{i}", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5),
Hist("ThetaVsP_d1", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz") 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, config, 3, 2, "daug1")
plotn(f, args, config, 3, 2, f"daug{i}")
config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, r = 5),
Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4), config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, xlo = 0, xhi = 2000),
Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.05), Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4, r = 10),
Hist("Pt_M", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 5), Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 0.2, ndiv = 5),
Hist("Phi_M", xtitle = "#phi [rad]"), Hist("Pt_M", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 50),
Hist("ThetaVsP_M", xtitle = "p^{0} [GeV]", ytitle = "#theta [rad]", logx = True, logy = True, d = "colz") 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"), # config = [Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
Hist("ThetaVsP_d1", xtitle = "p^{0} [GeV]", ytitle = "#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, config, 2, 2, "twod") #
# plotn(f, args, config, 2, 2, "twod")
config = [Hist("Vtx_X", xtitle = "x [mm]", r = 5), config = [Hist("Vtx_X", xtitle = "x [mm]", r = 5),
Hist("Vtx_Y", xtitle = "y [mm]", r = 5), Hist("Vtx_Y", xtitle = "y [mm]", r = 5),
Hist("Vtx_Z", xtitle = "z [mm]", r = 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_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")
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment