From 84a463298706a19cef2e0a1eb8268bc13a2373b5 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Mon, 19 Dec 2022 10:51:16 +0000
Subject: [PATCH] Update foresee sample config to 13.6 TeV and add ability to
sample hists in 1 or 2D
---
.../Generation/python/faser_parser.py | 14 ++++++--
.../share/generate_forsee_events.py | 4 +--
.../ForeseeGenerator/share/plot_validation.py | 1 +
.../ForeseeGenerator/share/validate_grid.py | 6 ++--
Generators/ParticleGun/python/histsampling.py | 16 ++++-----
Generators/ParticleGun/python/samplers.py | 35 +++++++++++++++----
6 files changed, 54 insertions(+), 22 deletions(-)
diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 43f98387e..ff8f1d379 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -4,6 +4,13 @@
#
# Parser function for particle gun samples
#
+
+def float_or_none(arg):
+ if arg is None or arg == "None":
+ return None
+
+ return float(arg)
+
def faser_pgparser():
import sys
@@ -35,13 +42,16 @@ def faser_pgparser():
help="Specify particle mass (in MeV)")
parser.add_argument("--radius", default=100., type=float,
help="Specify radius (in mm)")
- parser.add_argument("--angle", default=0.005, type=float,
+ parser.add_argument("--angle", default=0.005, type=float_or_none,
help="Specify angular width (in Rad)")
parser.add_argument("--zpos", default=None, type=float,
help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
parser.add_argument("--sampler", default="log",
- help="Specify energy sampling (log, lin, const)")
+ help="Specify energy sampling (log, lin, const, hist, hist2D)")
+ parser.add_argument("--hist_name", default="log",
+ help="Specify energy sampling name for hist sampler file.root:hist")
+
parser.add_argument("--minE", default=10., type=float,
help="Minimum energy in GeV (for log or lin sampler)")
parser.add_argument("--maxE", default=1000., type=float,
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 3d84bb153..0c3096f8a 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -459,7 +459,7 @@ def production(args):
f.write()
if args.outdir:
- p = proc.run(["cp", filename, outdir + "/" + filename.split("/")[-1]])
+ p = proc.run(["cp", filename, outdir + "/" + filename.split("/")[-1].replace(".hepmc", f"{nrun}.hepmc")])
if not p.returncode:
proc.run(["rm", filename])
@@ -492,7 +492,7 @@ if __name__ == "__main__":
parser.add_argument("--couplings", "-c", required = False, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
parser.add_argument("--pid1", default = None, required = False, type = int, help = "PID of daughter 1")
parser.add_argument("--pid2", default = None, type = int, help = "PID of daughter 2 (if not set then will be -PID1)")
- parser.add_argument("--Ecom", default = "14", help = "Center of mass energy [TeV]")
+ parser.add_argument("--Ecom", default = "13.6", help = "Center of mass energy [TeV]")
parser.add_argument("--outdir", "-o", default = None, help = "Output path")
parser.add_argument("--path", default = ".", help = "Path to foresee installation")
parser.add_argument("--hepmc", action = "store_true", help = "Write HepMC events")
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index d4ff9c8b2..b3cfc2b55 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -65,6 +65,7 @@ def plotn(f, args, configs, x, y, outname):
c._objs.append(plot(f, *cfg))
c.Print(f"{args.output}/{outname}.eps")
+ c.Print(f"{args.output}/{outname}.png")
return
diff --git a/Generators/ForeseeGenerator/share/validate_grid.py b/Generators/ForeseeGenerator/share/validate_grid.py
index e3d5302df..cf0146848 100644
--- a/Generators/ForeseeGenerator/share/validate_grid.py
+++ b/Generators/ForeseeGenerator/share/validate_grid.py
@@ -7,15 +7,15 @@ grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
with open(grid) as f:
data = json.load(f)
-name = "DarkPhotonGrid2"
+name = "DarkPhotonGrid13p6"
path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
-ecom = 14
+ecom = 13.6
pid1 = 11
pid2 = -11
for coup, masses in data["samples"].items():
for m in masses:
- infile = f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):e}to_{pid1}_{pid2}_{name}.hepmc"
+ infile = f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}.hepmc"
outfile = infile.replace(".hepmc", ".EVNT.pool.root")
valfile = infile.replace(".hepmc", ".VAL.pool.root")
valdir = infile.replace(".hepmc", "")
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
index c64112cc8..9f541da7c 100644
--- a/Generators/ParticleGun/python/histsampling.py
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -27,6 +27,9 @@ def load_hist(*args):
h = args[0].Get(args[1]).Clone()
if h is None:
raise Exception("Error in histogram loading from " + args)
+ else:
+ h.SetDirectory(0)
+
return h
@@ -40,7 +43,10 @@ def get_sampling_vars(h):
globalbin_to_axisbin = {} # for reverse axis bin lookup to get edges
globalbins = [] # because they aren't easily predicted, nor contiguous
cheights = [0] # cumulative "histogram" from which to uniformly sample
- if issubclass(type(h), ROOT.TH1):
+
+
+
+ if issubclass(type(h), ROOT.TH1) and not issubclass(type(h), ROOT.TH2):
for ix in range(1, h.GetNbinsX()+1):
iglobal = h.GetBin(ix)
globalbins.append(iglobal)
@@ -109,10 +115,6 @@ class TH1(object):
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th1)
return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin)
- def __getattr__(self, attr):
- "Forward all attributes to the contained TH1"
- return getattr(self.th1, attr)
-
class TH2(object):
"Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
@@ -127,6 +129,4 @@ class TH2(object):
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th2)
return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin)
- def __getattr__(self, attr):
- "Forward other attributes to the contained TH2"
- return getattr(self.th2, attr)
+
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index 90e9676a5..d6e540ea6 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -1,7 +1,7 @@
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
import ROOT, math, random
-from ParticleGun.histsampling import TH1
+from ParticleGun.histsampling import TH1,TH2
## For convenience
PI = math.pi
@@ -191,13 +191,26 @@ class TH1Sampler(ContinuousSampler):
def __init__(self, *args):
self.hist = TH1(*args)
- if self.hist.GetEntries() < 1:
- raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.GetName())
+ if self.hist.th1.GetEntries() < 1:
+ raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.th1.GetName())
def shoot(self):
return self.hist.GetRandom()
+class TH2Sampler(ContinuousSampler):
+ "Randomly sample from a 2D ROOT histogram."
+
+ def __init__(self, *args):
+ self.hist = TH2(*args)
+ if self.hist.th2.GetEntries() < 1:
+ raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.th2.GetName())
+
+ def shoot(self):
+ return self.hist.GetRandom()
+
+
+
########################################
@@ -572,9 +585,12 @@ class EThetaMPhiSampler(MomSampler):
# TODO: ensure that E >= m!
- def __init__(self, energy, theta, mass=0.0, phi=[0, TWOPI]):
+ def __init__(self, energy, theta=None, mass=0.0, phi=[0, TWOPI]):
self.energy = energy
- self.theta = theta
+ if theta is None:
+ self._theta = None
+ else:
+ self.theta = theta
self.mass = mass
self.phi = phi
@@ -616,10 +632,15 @@ class EThetaMPhiSampler(MomSampler):
pz = p cos(theta)
pt = p sin(theta)
"""
- e = self.energy()
+
+ if self._theta is None:
+ e,theta = self.energy()
+ else:
+ e = self.energy()
+ theta = self.theta()
+
m = self.mass()
p = math.sqrt( e**2 - m**2 )
- theta = self.theta()
pz = p * math.cos(theta)
pt = p * math.sin(theta)
phi = self.phi()
--
GitLab