diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index 814235e363a8ef5fb03df7a266676ff80247c672..dd599bf846d4d00380a045ed9a543123d48c8405 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -34,6 +34,7 @@ class RadialPosSampler(Sampler):
@property
def r(self):
"r position sampler"
+
fwhm = 2*self.radius
sig = fwhm/(2 * sqrt(2 * log(2)))
@@ -44,8 +45,6 @@ class RadialPosSampler(Sampler):
y = random.gauss(0, self.radius)
return sqrt(x**2 + y**2)
- # return random.gauss(0, sig)
-
@property
def phi(self):
"phi position sampler"
@@ -62,7 +61,7 @@ class RadialPosSampler(Sampler):
return ROOT.TLorentzVector(x, y, z, t)
if __name__ == "__main__":
-# test when run from command line
+ # test when run from command line
import numpy as np
import matplotlib.pyplot as plt
@@ -70,8 +69,8 @@ if __name__ == "__main__":
xlist = []
ylist = []
- r = RadialPosSampler(r = 1, x = 10, y = -10, z = 10)
- for i in range (10000):
+ r = RadialPosSampler(r = -10, x = 0, y = 0, z = 10)
+ for i in range (100000):
res = r.shoot()
xlist.append(res.X())
ylist.append(res.Y())
@@ -81,10 +80,13 @@ if __name__ == "__main__":
plt.figure(figsize = (5,5))
plt.subplot(2,2,1)
- plt.scatter(xarr, yarr, marker = "*")
+ plt.hist2d(xarr, yarr)
plt.subplot(2,2,2)
plt.hist(yarr)
plt.subplot(2,2,3)
plt.hist(xarr)
+ plt.subplot(2,2,4)
+ plt.hist(np.sqrt(xarr**2 + yarr**2))
plt.tight_layout()
plt.show()
+x
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 9c6be63475b6e066ece1ceac0aea059b8837a0dc..0112555e0d384f7d37ee8e2b93cc7bb074775fe2 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"
@@ -80,7 +81,7 @@ class EvgenValidation(EvgenAnalysisAlg):
def initialize(self):
# All daughters
- self.hists.add("PIDs", 60, -30, 30)
+ self.hists.add("PIDs", 600, -300, 300)
# Daughter i
tbins, pbins = self.binning()
@@ -92,26 +93,30 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists.add(f"Theta_d{i}", 20, 0, 0.001)
self.hists.add(f"Phi_d{i}", 16, -3.2, 3.2)
self.hists.add(f"ThetaVsP_d{i}", arrayX = tbins, arrayY = pbins)
- self.hists.add(f"Mass_d{i}", 200, 0, 0.01)
+ self.hists.add(f"Mass_d{i}", 5000, 0, 1)
# 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_XY", 50, -100, 100, 50, -100, 100)
+ self.hists.add("Vtx_X_LLP", 50, -100, 100)
+ self.hists.add("Vtx_Y_LLP", 50, -100, 100)
+ self.hists.add("Vtx_Z_LLP", 500, -5000, 0)
+ self.hists.add("Vtx_R_LLP", 20, 0, 200)
+ self.hists.add("Vtx_XY_LLP", 50, -100, 100, 50, -100, 100)
+
+ self.hists.add("Vtx_X_All", 50, -100, 100)
+ self.hists.add("Vtx_Y_All", 50, -100, 100)
+ self.hists.add("Vtx_Z_All", 500, -5000, 0)
+ self.hists.add("Vtx_R_All", 20, 0, 200)
+ self.hists.add("Vtx_XY_All", 50, -100, 100, 50, -100, 100)
return StatusCode.Success
@@ -137,11 +142,12 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists["PIDs"].Fill(p.pdg_id(), self.weight)
return
- def fillVertex(self, v):
- 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)
+ def fillVertex(self, label, v):
+ self.hists[f"Vtx_X_{label}"].Fill(v.x(), self.weight)
+ self.hists[f"Vtx_Y_{label}"].Fill(v.y(), self.weight)
+ self.hists[f"Vtx_Z_{label}"].Fill(v.z(), self.weight)
+ self.hists[f"Vtx_XY_{label}"].Fill(v.x(), v.y(), self.weight)
+ self.hists[f"Vtx_R_{label}"].Fill(sqrt(v.x()**2 + v.y()**2), self.weight)
return
@@ -151,26 +157,38 @@ class EvgenValidation(EvgenAnalysisAlg):
# Loop over all particles in events (assuming mother not stored)
momenta = []
+ vertices = []
mother = HepMC.FourVector(0,0,0,0)
llp_vtx = None
for i, p in enumerate(evt.particles):
#p.print()
self.fillDaughter(p)
+
momenta.append(p.momentum())
mother += p.momentum()
- if i == 0:
+
+ if i == 0 and p.production_vertex():
#p.production_vertex().print()
llp_vtx = p.production_vertex().point3d()
+ if p.production_vertex():
+ vertices.append(p.production_vertex().point3d())
+
# Fill daughter plots
for i in range(self.ndaughters):
+ if i >= len(momenta): continue
self.fillKin(f"d{i}", momenta[i])
# Fill mother plots
self.fillKin("M", mother, mass = True)
- # Fill vertex plots
- self.fillVertex(llp_vtx)
+ # Fill all vertices
+ for v in vertices:
+ self.fillVertex("All", v)
+
+ # Fill LLP vertex plots
+ if llp_vtx:
+ self.fillVertex("LLP", llp_vtx)
return StatusCode.Success
@@ -184,7 +202,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 +236,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..49436e95fe348844da228d7eff38059139a37b40 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,88 @@ 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")
+ 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.1, 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 = 10000),
+ Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4, r = 10),
+ Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 1., 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, "daug0")
+ plotn(f, args, config, 3, 2, "mother")
- 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, args, Hist("PIDs", xtitle="PDG Id", xlo=-30, xhi=30), 1, 1, "pid")
+
+
+# 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_LLP", xtitle = "x [mm]", r = 5),
+ Hist("Vtx_Y_LLP", xtitle = "y [mm]", r = 5),
+ Hist("Vtx_Z_LLP", xtitle = "z [mm]", r = 5, ndiv = 5),
+ Hist("Vtx_XY_LLP", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
+ Hist("Vtx_R_LLP", xtitle = "r [mm]", r = 5, ndiv = 5)
]
- plotn(f, config, 3, 2, "daug1")
+ plotn(f, args, config, 3, 2, "vtx_llp")
- 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")
- ]
- plotn(f, config, 3, 2, "mother")
+ config = [Hist("Vtx_X_All", xtitle = "x [mm]", r = 5),
+ Hist("Vtx_Y_All", xtitle = "y [mm]", r = 5),
+ Hist("Vtx_Z_All", xtitle = "z [mm]", r = 5, ndiv = 5),
+ Hist("Vtx_XY_All", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
+ Hist("Vtx_R_All", xtitle = "r [mm]", r = 5, ndiv = 5)
+ ]
- plotn(f, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid")
+ plotn(f, args, config, 3, 2, "vtx_all")
- 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("Vtx_X", xtitle = "x [mm]", r = 5),
+# Hist("Vtx_Y", xtitle = "y [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_R", xtitle = "r [mm]", r = 5, ndiv = 5)
+# ]
+#
+# plotn(f, args, config, 3, 2, "vtx")
- 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))
- ]
- plotn(f, config, 2, 2, "vtx")
+
diff --git a/Generators/GeneratorUtils/python/ShiftLOSConfig.py b/Generators/GeneratorUtils/python/ShiftLOSConfig.py
index 4d7f02d59d93591aa6c0e26156ce9a07e10e445d..b9c3ec2a5a209e9c555739099128a44d691bf87e 100644
--- a/Generators/GeneratorUtils/python/ShiftLOSConfig.py
+++ b/Generators/GeneratorUtils/python/ShiftLOSConfig.py
@@ -7,6 +7,8 @@ from AthenaConfiguration.MainServicesConfig import AthSequencer
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator, ConfigurationError
from AthenaConfiguration.ComponentFactory import CompFactory
+from SGComps.AddressRemappingConfig import InputOverwriteCfg
+
from GeneratorUtils.ShiftLOS import ShiftLOS
@@ -14,9 +16,13 @@ def ShiftLOSCfg(ConfigFlags, **kwargs) :
import McParticleEvent.Pythonizations
cfg = ComponentAccumulator()
+
+ # Rename old truth collection to add ATLAS coord to can still use BeamTruthEvent for the one in FASER Coords
+ cfg.merge(InputOverwriteCfg("McEventCollection", "BeamTruthEvent", "McEventCollection", "BeamTruthEvent_ATLASCoord"))
+
shift = ShiftLOS(name = kwargs.setdefault("name", "ShiftLOS"))
- shift.InputMCEventKey = kwargs.setdefault("InputMCEventKey", "BeamTruthEvent")
- shift.OutputMCEventKey = kwargs.setdefault("OutputMCEventKey", "BeamTruthEventShifted")
+ shift.InputMCEventKey = kwargs.setdefault("InputMCEventKey", "BeamTruthEvent_ATLASCoord")
+ shift.OutputMCEventKey = kwargs.setdefault("OutputMCEventKey", "BeamTruthEvent")
shift.xcross = kwargs.setdefault("xcross", 0)
shift.ycross = kwargs.setdefault("ycross", 0)
shift.xshift = kwargs.setdefault("xshift", 0)
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index 156f20d22770b1aef4a0daf8ea117833e1b2cdef..6b99fe63410583ff57260d48b813ba285949d906 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -137,21 +137,20 @@ if __name__ == '__main__':
if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift):
-
- MCEventKey = "BeamTruthEventShifted"
+
import McParticleEvent.Pythonizations
from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg
- cfg.merge(ShiftLOSCfg(ConfigFlags, OutputMCEventKey = MCEventKey,
+ cfg.merge(ShiftLOSCfg(ConfigFlags,
xcross = ConfigFlags.Sim.Beam.xangle, ycross = ConfigFlags.Sim.Beam.yangle,
xshift = ConfigFlags.Sim.Beam.xshift, yshift = ConfigFlags.Sim.Beam.yshift))
- else:
- MCEventKey = "BeamTruthEvent"
#
# Add the G4FaserAlg
#
from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg
- cfg.merge(G4FaserAlgCfg(ConfigFlags, InputTruthCollection = MCEventKey))
+ cfg.merge(G4FaserAlgCfg(ConfigFlags))
+
+ ###cfg.getEventAlgo("OutputStreamHITS").ItemList += ["McEventCollection#BeamTruthEvent_ATLASCoord"]
#
# Dump config
#