From 51c5ec1d7d1cfa2d20a62c6da3dcd43b03ee0c5e Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Thu, 26 May 2022 18:39:31 +0100
Subject: [PATCH 01/11] Update validation scripts
---
.../ForeseeGenerator/python/Validate.py | 26 +++---
.../ForeseeGenerator/share/plot_validation.py | 90 ++++++++++---------
2 files changed, 61 insertions(+), 55 deletions(-)
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 9c6be6347..8f27b36ab 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 8d5a04701..5f702f695 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")
--
GitLab
From 9da3fc6a2944722dd7aac3da97997d5283b52739 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Thu, 26 May 2022 19:05:23 +0100
Subject: [PATCH 02/11] Bug fix for radial sampler
---
.../python/RadialPosSampler.py | 20 ++++++++++++-------
1 file changed, 13 insertions(+), 7 deletions(-)
diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index a0ae101c5..ad7d4f273 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -34,14 +34,18 @@ class RadialPosSampler(Sampler):
@property
def r(self):
"r position sampler"
+
fwhm = 2*self.radius
sig = fwhm/(2 * sqrt(2 * log(2)))
if self.radius < 0:
- return random.uniform(0, abs(self.radius))
+ return sqrt(random.uniform(0, abs(self.radius**2)))
else:
- return random.gauss(0, self.radius)
- # return random.gauss(0, sig)
+ x = random.gauss(0, self.radius)
+ y = random.gauss(0, self.radius)
+ return sqrt(x**2 + y**2)
+
+ #return random.gauss(0, sig)
@property
def phi(self):
@@ -59,7 +63,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
@@ -67,8 +71,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())
@@ -78,10 +82,12 @@ 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()
--
GitLab
From 10c85b069f560abe5323b3c07197afb69f75a936 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 27 May 2022 18:56:38 +0100
Subject: [PATCH 03/11] Deal with renaming the BeamTruthEvent in ShiftLOS
gracefully
---
Generators/GeneratorUtils/python/ShiftLOSConfig.py | 10 ++++++++--
.../G4FaserAlg/test/G4FaserAlgConfigNew_Test.py | 11 +++++------
2 files changed, 13 insertions(+), 8 deletions(-)
diff --git a/Generators/GeneratorUtils/python/ShiftLOSConfig.py b/Generators/GeneratorUtils/python/ShiftLOSConfig.py
index 4d7f02d59..b9c3ec2a5 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 156f20d22..6b99fe634 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
#
--
GitLab
From a40f55b918b946d9023e102ca6761d2d1a0c6977 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Tue, 31 May 2022 15:04:45 +0100
Subject: [PATCH 04/11] Update validation code
---
.../ForeseeGenerator/python/Validate.py | 50 +++++++++++++------
.../ForeseeGenerator/share/plot_validation.py | 43 ++++++++++++----
2 files changed, 67 insertions(+), 26 deletions(-)
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 8f27b36ab..0112555e0 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -81,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()
@@ -93,7 +93,7 @@ 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", 1000, 0, 10000)
@@ -106,11 +106,17 @@ class EvgenValidation(EvgenAnalysisAlg):
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)
- 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)
+ 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
@@ -136,12 +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)
- self.hists["Vtx_R"].Fill(sqrt(v.x()**2 + v.y()**2), 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
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index 5f702f695..49436e95f 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -90,7 +90,7 @@ if __name__ == "__main__":
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"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")
@@ -98,9 +98,9 @@ if __name__ == "__main__":
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),
+ 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 = 0.2, ndiv = 5),
+ 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")
@@ -108,7 +108,7 @@ if __name__ == "__main__":
plotn(f, args, config, 3, 2, "mother")
- plotn(f, args, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid")
+ 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"),
@@ -118,11 +118,34 @@ if __name__ == "__main__":
#
# 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, 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)
+ 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, args, config, 3, 2, "vtx")
+ plotn(f, args, config, 3, 2, "vtx_llp")
+
+
+ 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, args, config, 3, 2, "vtx_all")
+
+
+# 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")
+
+
+
--
GitLab
From 6a36c53a19e4bb5c906d63d3cba636db5697f37a Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 3 Jun 2022 10:27:33 +0100
Subject: [PATCH 05/11] Fix LOS dealing with ATLAS vs FASER coords
---
.../GeneratorUtils/python/ShiftLOSConfig.py | 5 ---
.../test/G4FaserAlgConfigNew_Test.py | 40 ++++++++++++++-----
2 files changed, 30 insertions(+), 15 deletions(-)
diff --git a/Generators/GeneratorUtils/python/ShiftLOSConfig.py b/Generators/GeneratorUtils/python/ShiftLOSConfig.py
index b9c3ec2a5..ffc5c79b9 100644
--- a/Generators/GeneratorUtils/python/ShiftLOSConfig.py
+++ b/Generators/GeneratorUtils/python/ShiftLOSConfig.py
@@ -7,8 +7,6 @@ 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
@@ -17,9 +15,6 @@ def ShiftLOSCfg(ConfigFlags, **kwargs) :
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_ATLASCoord")
shift.OutputMCEventKey = kwargs.setdefault("OutputMCEventKey", "BeamTruthEvent")
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index 6b99fe634..f56b7bf73 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -58,11 +58,21 @@ if __name__ == '__main__':
import sys
ConfigFlags.fillFromArgs(sys.argv[1:])
-# from math import atan
-# from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m
-# from AthenaCommon.PhysicalConstants import pi
-# import ParticleGun as PG
-# ConfigFlags.Sim.Gun = {"Generator" : "SingleParticle", "pid" : 11, "energy" : PG.LogSampler(10*GeV, 1*TeV), "theta" : PG.GaussianSampler(0, atan((10*cm)/(7*m)), oneside = True), "phi" : [0, 2*pi], "mass" : 0.511, "radius" : -10*cm, "randomSeed" : 12345}
+ doShiftLOS = (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
+ ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift)
+
+ from math import atan
+ from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m
+ from AthenaCommon.PhysicalConstants import pi
+ import ParticleGun as PG
+ ConfigFlags.Sim.Gun = {"Generator" : "SingleParticle", "pid" : 11, "energy" : PG.LogSampler(10*GeV, 1*TeV), "theta" :
+ PG.GaussianSampler(0, atan((10*cm)/(7*m)), oneside = True), "phi" : [0, 2*pi], "mass" : 0.511, "radius" : -10*cm, "randomSeed" : 12345}
+
+ if doShiftLOS:
+ pgConfig = ConfigFlags.Sim.Gun
+ pgConfig["McEventKey"] = "BeamTruthEvent_ATLASCoord"
+ ConfigFlags.Sim.Gun = pgConfig
+
#
# By being a little clever, we can steer the geometry setup from the command line using GeoModel.FaserVersion
@@ -97,12 +107,16 @@ if __name__ == '__main__':
print("Input.Files = ",ConfigFlags.Input.Files)
#
-# If so, and only one file that ends in .events read as HepMC
+# If so, and only one file that ends in .events or .hepmc read as HepMC
#
if len(ConfigFlags.Input.Files) == 1 and (ConfigFlags.Input.Files[0].endswith(".events") or ConfigFlags.Input.Files[0].endswith(".hepmc")):
from HEPMCReader.HepMCReaderConfig import HepMCReaderCfg
- cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ if doShiftLOS:
+ cfg.merge(HepMCReaderCfg(ConfigFlags, McEventKey = "BeamTruthEvent_ATLASCoord"))
+ else:
+ cfg.merge(HepMCReaderCfg(ConfigFlags))
from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
cfg.merge(McEventSelectorCfg(ConfigFlags))
@@ -113,6 +127,12 @@ if __name__ == '__main__':
else:
from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
cfg.merge(PoolReadCfg(ConfigFlags))
+
+ if doShiftLOS:
+ from SGComps.AddressRemappingConfig import InputOverwriteCfg
+ # 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"))
+
#
# If not, configure the particle gun as requested, or using defaults
#
@@ -122,6 +142,7 @@ if __name__ == '__main__':
#
from FaserParticleGun.FaserParticleGunConfig import FaserParticleGunCfg
cfg.merge(FaserParticleGunCfg(ConfigFlags))
+
from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
cfg.merge(McEventSelectorCfg(ConfigFlags))
@@ -135,8 +156,7 @@ if __name__ == '__main__':
# Shift LOS
#
- if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
- ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift):
+ if doShiftLOS:
import McParticleEvent.Pythonizations
from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg
@@ -150,7 +170,7 @@ if __name__ == '__main__':
from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg
cfg.merge(G4FaserAlgCfg(ConfigFlags))
- ###cfg.getEventAlgo("OutputStreamHITS").ItemList += ["McEventCollection#BeamTruthEvent_ATLASCoord"]
+ #cfg.getEventAlgo("OutputStreamHITS").ItemList += ["McEventCollection#BeamTruthEvent_ATLASCoord"]
#
# Dump config
#
--
GitLab
From 7f77e7523b092e072ff91fb79ab61681b4c198fe Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 3 Jun 2022 10:28:04 +0100
Subject: [PATCH 06/11] Fix LOS dealing with ATLAS vs FASER coords
---
.../G4FaserAlg/test/G4FaserAlgConfigNew_Test.py | 12 ++++++------
1 file changed, 6 insertions(+), 6 deletions(-)
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index f56b7bf73..998b40c28 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -61,12 +61,12 @@ if __name__ == '__main__':
doShiftLOS = (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift)
- from math import atan
- from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m
- from AthenaCommon.PhysicalConstants import pi
- import ParticleGun as PG
- ConfigFlags.Sim.Gun = {"Generator" : "SingleParticle", "pid" : 11, "energy" : PG.LogSampler(10*GeV, 1*TeV), "theta" :
- PG.GaussianSampler(0, atan((10*cm)/(7*m)), oneside = True), "phi" : [0, 2*pi], "mass" : 0.511, "radius" : -10*cm, "randomSeed" : 12345}
+# from math import atan
+# from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m
+# from AthenaCommon.PhysicalConstants import pi
+# import ParticleGun as PG
+# ConfigFlags.Sim.Gun = {"Generator" : "SingleParticle", "pid" : 11, "energy" : PG.LogSampler(10*GeV, 1*TeV), "theta" :
+# PG.GaussianSampler(0, atan((10*cm)/(7*m)), oneside = True), "phi" : [0, 2*pi], "mass" : 0.511, "radius" : -10*cm, "randomSeed" : 12345}
if doShiftLOS:
pgConfig = ConfigFlags.Sim.Gun
--
GitLab
From 588fbca893687b8716819dbcc99412780794215d Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 3 Jun 2022 12:18:10 +0100
Subject: [PATCH 07/11] Fix skipping events with HepMC and add hepmc to evnt
converter
---
.../share/convert_hepmc_to_evnt.py | 62 +++++++++++++++++++
.../HEPMCReader/python/HepMCReaderConfig.py | 46 ++++++++++++--
2 files changed, 104 insertions(+), 4 deletions(-)
create mode 100644 Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
diff --git a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
new file mode 100644
index 000000000..88bda6736
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
@@ -0,0 +1,62 @@
+def getNEvents(fname, maxEvents):
+ "Work out how many events are in the file"
+
+ n = 0
+ with open(fname) as f:
+ n = sum(1 for l in f if l.startswith("E "))
+
+ if maxEvents != -1 and n > maxEvents:
+ n = maxEvents
+
+ print ("Setting Maximum events to", n)
+ return n
+
+
+if __name__ == "__main__":
+
+ import sys
+
+ from AthenaCommon.Logging import log
+ from AthenaCommon.Constants import DEBUG
+ log.setLevel(DEBUG)
+
+ from AthenaCommon.Configurable import Configurable
+ Configurable.configurableRun3Behavior = 1
+
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ ConfigFlags.Input.RunNumber = [12345]
+ ConfigFlags.Input.OverrideRunNumber = True
+ ConfigFlags.Input.LumiBlockNumber = [1]
+
+ ConfigFlags.Input.isMC = True
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion
+ ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Default FASER geometry
+ ConfigFlags.Detector.EnableFaserSCT = True
+
+ ConfigFlags.Output.EVNTFileName = "myEVNT.pool.root"
+
+ ConfigFlags.Exec.MaxEvents= -1
+
+ import sys
+ ConfigFlags.fillFromArgs(sys.argv[1:])
+
+
+ ConfigFlags.Exec.MaxEvents = getNEvents(ConfigFlags.Input.Files[0], ConfigFlags.Exec.MaxEvents)
+
+ ConfigFlags.lock()
+
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ cfg = MainServicesCfg(ConfigFlags)
+
+ from HEPMCReader.HepMCReaderConfig import HepMCReaderCfg
+ cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
+ cfg.merge(McEventSelectorCfg(ConfigFlags))
+
+ itemList = [ "McEventCollection#*" ]
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ cfg.merge(OutputStreamCfg(ConfigFlags, "EVNT", itemList, disableEventTag = True))
+
+ sc = cfg.run()
+ sys.exit(not sc.isSuccess())
diff --git a/Generators/HEPMCReader/python/HepMCReaderConfig.py b/Generators/HEPMCReader/python/HepMCReaderConfig.py
index 72081a86f..bb116c11c 100644
--- a/Generators/HEPMCReader/python/HepMCReaderConfig.py
+++ b/Generators/HEPMCReader/python/HepMCReaderConfig.py
@@ -2,21 +2,59 @@
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-# import sys
+import sys, tempfile
from AthenaConfiguration.MainServicesConfig import AthSequencer
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
from TruthIO.TruthIOConf import HepMCReadFromFile
+def get_file_skip_events(ConfigFlags):
+ "Create a file with events from ConfigFlags.Exec.SkipEvents to ConfigFlags.Exec.SkipEvents + ConfigFlags.Exec.MaxEvents"
+
+ usetemp = True
+
+ skip = ConfigFlags.Exec.SkipEvents
+ fname = ConfigFlags.Input.Files[0]
+ evtMax = ConfigFlags.Exec.MaxEvents
+
+ if skip == 0:
+ return fname
+
+ if usetemp:
+ fout = tempfile.NamedTemporaryFile("w", delete = False)
+ foutname = fout.name
+ else:
+ foutname, fext = ".".join(fname.split('.')[:-1]), fname.split('.')[-1]
+ foutname = f"{foutname}-evts{skip}-{skip+evtMax}.{fext}"
+ fout = open(foutname, "w")
+
+ fout.write("HepMC::Version 2.06.09\nHepMC::IO_GenEvent-START_EVENT_LISTING\n")
+
+ ievt = 0
+ with open(fname) as fin:
+ for l in fin:
+ if l.startswith("E "):
+ ievt += 1
+
+ if ievt > skip + evtMax:
+ break
+
+ if ievt > skip:
+ fout.write(l)
-def HepMCReaderCfg(ConfigFlags, **kwargs) :
- cfg = ComponentAccumulator(AthSequencer("AthBeginSeq", Sequential = True))
+ fout.write("HepMC::IO_GenEvent-END_EVENT_LISTING\n")
+ fout.close()
+
+ return foutname
+
+def HepMCReaderCfg(ConfigFlags, **kwargs) :
+ cfg = ComponentAccumulator(AthSequencer("AthBeginSeq", Sequential = True))
from TruthIO.TruthIOConf import HepMCReadFromFile
hepmc = CompFactory.HepMCReadFromFile(name = kwargs.setdefault("name", "FASERHepMCReader"))
- hepmc.InputFile = ConfigFlags.Input.Files[0]
+ hepmc.InputFile = get_file_skip_events(ConfigFlags) # ConfigFlags.Input.Files[0]
hepmc.McEventKey = kwargs.setdefault("McEventKey", "BeamTruthEvent")
cfg.addEventAlgo(hepmc, sequenceName = "AthBeginSeq", primary = True) # to run *before* G4
--
GitLab
From 3317fa17e584eff7e7b8e74937a16ed53ead5fc1 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 3 Jun 2022 12:19:21 +0100
Subject: [PATCH 08/11] Update validation code
---
.../ForeseeGenerator/python/Validate.py | 24 +++++++++++++------
.../ForeseeGenerator/share/plot_validation.py | 4 ++--
2 files changed, 19 insertions(+), 9 deletions(-)
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 0112555e0..399d85ae7 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -64,9 +64,10 @@ class HistSvc(object):
class EvgenValidation(EvgenAnalysisAlg):
"Gen-level validation"
- def __init__(self, name = "EvgenValidation", ndaughters = 2, outname = "validation.root"):
+ def __init__(self, name = "EvgenValidation", ndaughters = 2, outname = "validation.root", mother_stored = False):
super(EvgenValidation, self).__init__(name=name)
self.hists = HistSvc()
+ self.mother_stored = mother_stored
self.ndaughters = ndaughters
self.outname = outname
@@ -153,19 +154,27 @@ class EvgenValidation(EvgenAnalysisAlg):
def execute(self):
evt = self.events()[0]
- self.weight = evt.weights()[0]
+ self.weight = evt.weights()[0] if evt.weights() else 1
- # Loop over all particles in events (assuming mother not stored)
+ # Loop over all particles in events
momenta = []
vertices = []
mother = HepMC.FourVector(0,0,0,0)
llp_vtx = None
+
for i, p in enumerate(evt.particles):
- #p.print()
+ print("--- ", i)
+ p.print()
self.fillDaughter(p)
- momenta.append(p.momentum())
- mother += p.momentum()
+ if self.mother_stored:
+ if i == 0:
+ mother = p.momentum()
+ else:
+ momenta.append(p.momentum())
+ else:
+ momenta.append(p.momentum())
+ mother += p.momentum()
if i == 0 and p.production_vertex():
#p.production_vertex().print()
@@ -203,6 +212,7 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Run gen-level validation")
parser.add_argument("file", nargs="+", help = "full path to imput file")
parser.add_argument("--ndaughters", "-d", default = 2, type = int, help = "Number of daugthers to plot")
+ parser.add_argument("--mother_stored", "-m", default = False, action = "store_true", help = "Is mother stored in input?")
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")
@@ -236,7 +246,7 @@ if __name__ == "__main__":
fix()
acc = ComponentAccumulator()
- valid = EvgenValidation("EvgenValidation", ndaughters = args.ndaughters, outname = args.output)
+ valid = EvgenValidation("EvgenValidation", ndaughters = args.ndaughters, outname = args.output, mother_stored = args.mother_stored)
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 49436e95f..8429978bc 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -98,7 +98,7 @@ if __name__ == "__main__":
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),
+ config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, xlo = 0, xhi = 10000, r=10),
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),
@@ -108,7 +108,7 @@ if __name__ == "__main__":
plotn(f, args, config, 3, 2, "mother")
- plotn(f, args, Hist("PIDs", xtitle="PDG Id", xlo=-30, xhi=30), 1, 1, "pid")
+ plotn(f, args, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid")
# config = [Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
--
GitLab
From c5e058e7479a88dd64ff2fc70d718a7ca2d6f05d Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 3 Jun 2022 15:07:54 +0100
Subject: [PATCH 09/11] Fix skipping events with HepMC and add hepmc to evnt
converter
---
Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
index 88bda6736..382ad439e 100644
--- a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
+++ b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
@@ -54,7 +54,7 @@ if __name__ == "__main__":
from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
cfg.merge(McEventSelectorCfg(ConfigFlags))
- itemList = [ "McEventCollection#*" ]
+ itemList = [ "EventInfo#McEventInfo", "McEventCollection#*" ]
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
cfg.merge(OutputStreamCfg(ConfigFlags, "EVNT", itemList, disableEventTag = True))
--
GitLab
From cb51b400c503fb634c26ddf9086eb9935566cdee Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Tue, 7 Jun 2022 11:01:15 +0100
Subject: [PATCH 10/11] Merge
---
...serMC-MDC_PG_muon_fasernu_logE-101302.json | 2 +-
.../Generation/scripts/faserMDC_foresee.py | 20 ++++++-----
.../scripts/faserMDC_particlegun.py | 20 ++++++-----
.../Reconstruction/scripts/faserMDC_reco.py | 2 +-
.../scripts/submit_faserMDC_reco.sh | 2 +-
.../Simulation/scripts/faserMDC_simulate.py | 34 ++++++++++++-------
.../scripts/submit_faserMDC_simulate.sh | 30 ++++++++++++++--
.../python/RadialPosSampler.py | 1 -
.../HEPMCReader/python/HepMCReaderConfig.py | 30 +++++++++++-----
.../test/G4FaserAlgConfigNew_Test.py | 0
10 files changed, 98 insertions(+), 43 deletions(-)
mode change 100644 => 100755 Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101302.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101302.json
index 8ae524dc0..76ba19e11 100644
--- a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101302.json
+++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_muon_fasernu_logE-101302.json
@@ -9,5 +9,5 @@
"sampler": "log",
"segment": 0,
"short": "MDC_PG_muon_fasernu_logE",
- "zpos": -5000.0
+ "zpos": -4000.0
}
diff --git a/Control/CalypsoExample/Generation/scripts/faserMDC_foresee.py b/Control/CalypsoExample/Generation/scripts/faserMDC_foresee.py
index 607ce9046..b48bbf861 100755
--- a/Control/CalypsoExample/Generation/scripts/faserMDC_foresee.py
+++ b/Control/CalypsoExample/Generation/scripts/faserMDC_foresee.py
@@ -113,6 +113,14 @@ if __name__ == '__main__':
#if args.zpos:
# ConfigFlags.Sim.Gun["z"] = args.zpos
+ doShiftLOS = (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
+ ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift)
+
+ if doShiftLOS:
+ pgConfig = ConfigFlags.Sim.Gun
+ pgConfig["McEventKey"] = "BeamTruthEvent_ATLASCoord"
+ ConfigFlags.Sim.Gun = pgConfig
+
#
# MDC geometry configuration
#
@@ -153,26 +161,22 @@ if __name__ == '__main__':
# Shift LOS
#
- if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
- ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift):
-
- MCEventKey = "BeamTruthEventShifted"
+ if doShiftLOS:
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))
#
# Dump config
#
diff --git a/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py b/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py
index c9f9154ef..20970fdf3 100755
--- a/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py
@@ -124,6 +124,14 @@ if __name__ == '__main__':
# Pass this in one go to ConfigFlags
ConfigFlags.Sim.Gun = sg_dict
+ doShiftLOS = (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
+ ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift)
+
+ if doShiftLOS:
+ pgConfig = ConfigFlags.Sim.Gun
+ pgConfig["McEventKey"] = "BeamTruthEvent_ATLASCoord"
+ ConfigFlags.Sim.Gun = pgConfig
+
#
# MDC geometry configuration
#
@@ -164,26 +172,22 @@ if __name__ == '__main__':
# Shift LOS
#
- if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
- ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift):
-
- MCEventKey = "BeamTruthEventShifted"
+ if doShiftLOS:
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))
#
# Dump config
#
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
index f9003eece..942340f12 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
@@ -191,7 +191,7 @@ itemList = [ "xAOD::EventInfo#*"
#
if args.isMC:
# Add truth records here?
- itemList.extend( ["McEventCollection#*"] )
+ itemList.extend( ["McEventCollection#*", "TrackerSimDataCollection#*"] )
acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh
index 60c39ca37..b04d56aa9 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh
+++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh
@@ -124,7 +124,7 @@ fi
cd "$file_stem"
#
# Run job
-if [[ -z "$rtag" ]]; then
+if [[ -z "$tag" ]]; then
faserMDC_reco.py "--nevents=$nevents" "$file_path"
else
faserMDC_reco.py "--nevents=$nevents" "--reco=$tag" "$file_path"
diff --git a/Control/CalypsoExample/Simulation/scripts/faserMDC_simulate.py b/Control/CalypsoExample/Simulation/scripts/faserMDC_simulate.py
index 323adda69..12de492f4 100755
--- a/Control/CalypsoExample/Simulation/scripts/faserMDC_simulate.py
+++ b/Control/CalypsoExample/Simulation/scripts/faserMDC_simulate.py
@@ -127,6 +127,7 @@ if __name__ == '__main__':
# Skip events
#
ConfigFlags.Exec.SkipEvents = args.skip
+ ConfigFlags.Exec.MaxEvents = args.nevents
#
# Output file name
#
@@ -153,6 +154,9 @@ if __name__ == '__main__':
# import sys
# ConfigFlags.fillFromArgs(sys.argv[1:])
+ doShiftLOS = (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
+ ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift)
+
#
# MDC geometry configuration
#
@@ -178,7 +182,11 @@ if __name__ == '__main__':
if ConfigFlags.Input.Files[0].endswith(".events") or ConfigFlags.Input.Files[0].endswith(".hepmc"):
from HEPMCReader.HepMCReaderConfig import HepMCReaderCfg
- cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ if doShiftLOS:
+ cfg.merge(HepMCReaderCfg(ConfigFlags, McEventKey = "BeamTruthEvent_ATLASCoord"))
+ else:
+ cfg.merge(HepMCReaderCfg(ConfigFlags))
from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
cfg.merge(McEventSelectorCfg(ConfigFlags))
@@ -189,6 +197,11 @@ if __name__ == '__main__':
from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
cfg.merge(PoolReadCfg(ConfigFlags))
+ if doShiftLOS:
+ from SGComps.AddressRemappingConfig import InputOverwriteCfg
+ # 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"))
+
#
# Output file
#
@@ -198,26 +211,22 @@ if __name__ == '__main__':
#
# Shift LOS
#
- if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or
- ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift):
+ if doShiftLOS:
- 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))
#
# Dump config
#
@@ -235,10 +244,11 @@ if __name__ == '__main__':
# Execute and finish
#
- if args.verbose:
- cfg.foreach_component("*").OutputLevel = "DEBUG"
- else:
- cfg.foreach_component("*").OutputLevel = "INFO"
+ # This fails with ShiftLOSCfg...
+ #if args.verbose:
+ # cfg.foreach_component("*").OutputLevel = "DEBUG"
+ #else:
+ # cfg.foreach_component("*").OutputLevel = "INFO"
sc = cfg.run(maxEvents=args.nevents)
diff --git a/Control/CalypsoExample/Simulation/scripts/submit_faserMDC_simulate.sh b/Control/CalypsoExample/Simulation/scripts/submit_faserMDC_simulate.sh
index e00adfd50..3c175d725 100755
--- a/Control/CalypsoExample/Simulation/scripts/submit_faserMDC_simulate.sh
+++ b/Control/CalypsoExample/Simulation/scripts/submit_faserMDC_simulate.sh
@@ -2,7 +2,10 @@
# Used with a condor file to submit to vanilla universe
#
# Usage:
-# submit_faserMDC_simluate.sh input_file output_file [release_directory] [working_directory] [skip] [nevts]
+# submit_faserMDC_simluate.sh [--shift] input_file output_file [release_directory] [working_directory] [skip] [nevts]
+#
+# Options:
+# --shift - apply crossing angle (and FASER shift)
#
# input_file - full file name (with path)
# output_file - full output file name
@@ -22,6 +25,24 @@
SECONDS=0
#
# Parse command-line options
+while [ -n "$1" ]
+do
+ case "$1" in
+ -s | --shift)
+ echo "Applying crossing-angle shift"
+ xangle=1
+ shift;; # This 'eats' the argument
+
+ -*)
+ echo "Unknown option $1"
+ shift;;
+
+ *) break;; # Not an option, don't shift
+ esac
+done
+
+#
+# Parse command-line arguments
infile=${1}
outfile=${2}
release_directory=${3}
@@ -143,8 +164,11 @@ cd "${file_stem}"
# Run job
#if [[ -z "$tag" ]]; then
#fi
-faserMDC_simulate.py --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
-
+if [[ -z "$xangle" ]]; then
+ faserMDC_simulate.py --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+else
+ faserMDC_simulate.py --yangle -0.000150 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+fi
#
# Print out ending time
date
diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index dd599bf84..7f1ed1799 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -89,4 +89,3 @@ if __name__ == "__main__":
plt.hist(np.sqrt(xarr**2 + yarr**2))
plt.tight_layout()
plt.show()
-x
diff --git a/Generators/HEPMCReader/python/HepMCReaderConfig.py b/Generators/HEPMCReader/python/HepMCReaderConfig.py
index bb116c11c..affee04ea 100644
--- a/Generators/HEPMCReader/python/HepMCReaderConfig.py
+++ b/Generators/HEPMCReader/python/HepMCReaderConfig.py
@@ -2,7 +2,7 @@
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-import sys, tempfile
+import sys, tempfile, pathlib
from AthenaConfiguration.MainServicesConfig import AthSequencer
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
@@ -13,7 +13,8 @@ def get_file_skip_events(ConfigFlags):
"Create a file with events from ConfigFlags.Exec.SkipEvents to ConfigFlags.Exec.SkipEvents + ConfigFlags.Exec.MaxEvents"
usetemp = True
-
+ #usetemp = False
+
skip = ConfigFlags.Exec.SkipEvents
fname = ConfigFlags.Input.Files[0]
evtMax = ConfigFlags.Exec.MaxEvents
@@ -21,12 +22,22 @@ def get_file_skip_events(ConfigFlags):
if skip == 0:
return fname
+ print(f"get_file_skip_events skipping {skip} events with max {evtMax}")
+
if usetemp:
fout = tempfile.NamedTemporaryFile("w", delete = False)
foutname = fout.name
else:
- foutname, fext = ".".join(fname.split('.')[:-1]), fname.split('.')[-1]
- foutname = f"{foutname}-evts{skip}-{skip+evtMax}.{fext}"
+ infile = pathlib.Path(fname)
+ # Put this in current working directory
+ if evtMax > 0:
+ end = skip+evtMax
+ else:
+ end = 'all'
+
+ foutname = f"{infile.stem}=evts{skip}-{end}.{infile.suffix}"
+ #foutname, fext = ".".join(fname.split('.')[:-1]), fname.split('.')[-1]
+ #foutname = f"{foutname}-evts{skip}-{skip+evtMax}{fext}"
fout = open(foutname, "w")
fout.write("HepMC::Version 2.06.09\nHepMC::IO_GenEvent-START_EVENT_LISTING\n")
@@ -37,19 +48,22 @@ def get_file_skip_events(ConfigFlags):
if l.startswith("E "):
ievt += 1
- if ievt > skip + evtMax:
- break
+ if evtMax > 0 and ievt > skip + evtMax:
+ break
if ievt > skip:
+ #print(f"Writing event {ievt}")
fout.write(l)
-
+ # else:
+ # print(f"Skipping event {ievt}")
fout.write("HepMC::IO_GenEvent-END_EVENT_LISTING\n")
fout.close()
+ #print(f"Wrote to file {foutname}")
+
return foutname
-
def HepMCReaderCfg(ConfigFlags, **kwargs) :
cfg = ComponentAccumulator(AthSequencer("AthBeginSeq", Sequential = True))
from TruthIO.TruthIOConf import HepMCReadFromFile
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
old mode 100644
new mode 100755
--
GitLab
From 593f1b783fc38d5779116d230eb4570851b5c1f0 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Tue, 7 Jun 2022 11:04:27 +0100
Subject: [PATCH 11/11] Update foresee and fluka generation + validation
---
.../FlukaReader/python/FlukaReaderAlg.py | 20 +++---
.../ForeseeGenerator/python/Validate.py | 5 +-
.../share/generate_forsee_events.py | 61 +++++++++++++------
3 files changed, 55 insertions(+), 31 deletions(-)
diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py
index 3a8728e37..4cc449e65 100644
--- a/Generators/FlukaReader/python/FlukaReaderAlg.py
+++ b/Generators/FlukaReader/python/FlukaReaderAlg.py
@@ -333,8 +333,8 @@ class FlukaReader(EvgenAlg):
plt.figure()
ebins = np.linspace(0, 5000, 50)
plt.xlabel("Energy")
- plt.hist(self.before["E"], bins=ebins, histtype='step', color = "g", fill = False, label = "before")
- plt.hist(self.after["E"], bins = ebins, histtype='step', color = "r", fill = False, label = "after")
+ plt.hist(self.before["E"], bins=ebins, weights = np.array(self.before["w"])/self.nsamples, histtype='step', color = "g", fill = False, label = "before")
+ plt.hist(self.after["E"], bins = ebins, weights = np.array(self.after["w"])/self.nsamples, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("energy.png")
@@ -344,8 +344,8 @@ class FlukaReader(EvgenAlg):
thetaX = np.pi/2. - np.arccos(np.array(self.before["cosX"]))
thetaXout = np.pi/2. - np.arccos(np.array(self.after["cosX"]))
tbins = np.linspace(-0.5, 0.5, 100)
- plt.hist(thetaX, bins = tbins, histtype='step', color = "g", fill = False, label = "before")
- plt.hist(thetaXout, bins = tbins, histtype='step', color = "r", fill = False, label = "after")
+ plt.hist(thetaX, bins = tbins, weights = np.array(self.before["w"])/self.nsamples, histtype='step', color = "g", fill = False, label = "before")
+ plt.hist(thetaXout, bins = tbins, weights = np.array(self.after["w"])/self.nsamples, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("thetaX.png")
@@ -354,8 +354,8 @@ class FlukaReader(EvgenAlg):
plt.xlabel("Angle to beam in Y dir")
thetaY = np.pi/2. - np.arccos(np.array(self.before["cosY"]))
thetaYout = np.pi/2. - np.arccos(np.array(self.after["cosY"]))
- plt.hist(thetaY, bins = tbins, histtype='step', color = "g", fill = False, label = "before")
- plt.hist(thetaYout, bins = tbins, histtype='step', color = "r", fill = False, label = "after")
+ plt.hist(thetaY, bins = tbins, weights = np.array(self.before["w"])/self.nsamples, histtype='step', color = "g", fill = False, label = "before")
+ plt.hist(thetaYout, bins = tbins, weights = np.array(self.after["w"])/self.nsamples, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("thetaY.png")
@@ -363,16 +363,16 @@ class FlukaReader(EvgenAlg):
plt.figure()
plt.xlabel("Dispacement in X dir")
xbins = np.linspace(-300, 300, 100)
- plt.hist(self.before["x"], bins = xbins, histtype='step', color = "g", fill = False, label = "before")
- plt.hist(self.after["x"], bins = xbins, histtype='step', color = "r", fill = False, label = "after")
+ plt.hist(self.before["x"], bins = xbins, weights = np.array(self.before["w"])/self.nsamples, histtype='step', color = "g", fill = False, label = "before")
+ plt.hist(self.after["x"], bins = xbins, weights = np.array(self.after["w"])/self.nsamples, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("x.png")
plt.figure()
plt.xlabel("Dispacement in Y dir")
- plt.hist(self.before["y"], bins = xbins, histtype='step', color = "g", fill = False, label = "before")
- plt.hist(self.after["y"], bins = xbins, histtype='step', color = "r", fill = False, label = "after")
+ plt.hist(self.before["y"], bins = xbins, weights = np.array(self.before["w"])/self.nsamples, histtype='step', color = "g", fill = False, label = "before")
+ plt.hist(self.after["y"], bins = xbins, weights = np.array(self.after["w"])/self.nsamples, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("y.png")
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 399d85ae7..a93efeedd 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -163,8 +163,7 @@ class EvgenValidation(EvgenAnalysisAlg):
llp_vtx = None
for i, p in enumerate(evt.particles):
- print("--- ", i)
- p.print()
+ #p.print()
self.fillDaughter(p)
if self.mother_stored:
@@ -179,6 +178,8 @@ class EvgenValidation(EvgenAnalysisAlg):
if i == 0 and p.production_vertex():
#p.production_vertex().print()
llp_vtx = p.production_vertex().point3d()
+ elif i == 0 and p.end_vertex():
+ llp_vtx = p.end_vertex().point3d()
if p.production_vertex():
vertices.append(p.production_vertex().point3d())
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 37372b7bd..f02f4311f 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -19,7 +19,7 @@ class ForeseeGenerator(object):
self.daughter2_pid = daughter2_pid
self.outdir = outdir
self.path = path
- self.version = 1 # Forsee "version": 2 is in testing
+ self.version = 2 # Forsee "version": 2 is in testing
self.seed = randomSeed
# Set decay mode ...
@@ -102,6 +102,12 @@ class ForeseeGenerator(object):
coupling_ref=1,
massrange=[1.5, 10.]
)
+
+ self.model.set_br_1d(
+ modes = [self.mode],
+ filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"]
+ )
+
else:
masses_brem = [
0.01 , 0.0126, 0.0158, 0.02 , 0.0251, 0.0316, 0.0398,
@@ -134,6 +140,13 @@ class ForeseeGenerator(object):
masses = masses_dy,
)
+ self.model.set_br_1d(
+ modes = ["e_e", "mu_mu"],
+ finalstates = [[11, -11], [13, -13]],
+ filenames=["model/br/e_e.txt", "model/br/mu_mu.txt"],
+ #filenames=[f"model/br/all.txt"]
+ )
+
return self.decays()
@@ -175,6 +188,19 @@ class ForeseeGenerator(object):
nsample = 10,
)
+ if self.version == 1:
+ self.model.set_br_1d(
+ modes = [self.mode],
+ filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"]
+ )
+ else:
+ self.model.set_br_1d(
+ modes = ["gamma_gamma"],
+ finalstates = [[22, 22]],
+ filenames=["model/br/gamma_gamma.txt"]
+ #filenames=[f"model/br/all.txt"]
+ )
+
return self.decays()
@@ -185,24 +211,13 @@ class ForeseeGenerator(object):
self.model.set_ctau_1d(
filename=f"files/models/{self.modelname}/ctau.txt",
coupling_ref=1
- )
-
- self.model.set_br_1d(
- modes = [self.mode],
- filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"]
- )
+ )
else:
self.model.set_ctau_1d(
filename=f"model/ctau.txt",
coupling_ref=1
)
- self.model.set_br_1d(
- modes = [self.mode],
- finalstates = [[self.daughter1_pid, self.daughter2_pid]],
- filenames=[f"model/br/{self.mode}.txt"]
- )
-
# Get LLP spectrum
self.foresee.set_model(model=self.model)
# This is just a reference coupling
@@ -295,7 +310,15 @@ class ForeseeGenerator(object):
filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.hepmc"
- self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = -1500, seed = self.seed)
+ self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = -1.5, seed = self.seed, decaychannels = [self.mode])
+
+ cfgname = f"{self.foresee.dirpath}/Models/{self.modelname}/" + filename.replace(".hepmc", ".cfg")
+ print(f"save config to file: {cfgname}")
+ with open(cfgname, "w") as f:
+ f.write(" ".join(sys.argv))
+
+ return
+
def setup_foresee(path):
@@ -334,10 +357,10 @@ def add_to_python_path(path):
def parse_couplings(data, write_hepMC = False):
if write_hepMC:
- try:
- couplings = float(couplings)
- except ValueError:
- sus.exit("Only a single coupling allowed when writing HEPMC events")
+ if len(data) == 1:
+ couplings = float(data[0])
+ else:
+ sys.exit("Only a single coupling allowed when writing HEPMC events")
try:
couplings = [float(d) for d in data]
@@ -375,7 +398,7 @@ if __name__ == "__main__":
if args.pid2 is None:
args.pid2 = -args.pid1
- couplings = parse_couplings(args.couplings)
+ couplings = parse_couplings(args.couplings, args.hepmc)
print(f"Generating {args.model} events at Ecom = {args.Ecom}")
print(f" mother mass = {args.mass} GeV")
--
GitLab