From a81e5bd6ed6173e9b6e92a603840ff155a1dbf29 Mon Sep 17 00:00:00 2001
From: FaserMC <fasermc@cern.ch>
Date: Thu, 9 Mar 2023 18:21:47 +0100
Subject: [PATCH] Consolodate s0010 script updates
---
.../Generation/python/faser_parser.py | 4 ++
.../Generation/scripts/faser_particlegun.py | 45 ++++++++++++---
.../scripts/submit_faser_particlegun.sh | 39 +++++++++++--
Generators/ParticleGun/python/histsampling.py | 55 ++++++++++++++++---
Generators/ParticleGun/python/samplers.py | 23 +++++---
5 files changed, 137 insertions(+), 29 deletions(-)
diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 055dd7200..03817efb8 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -44,6 +44,10 @@ def faser_pgparser():
help="Specify radius (in mm)")
parser.add_argument("--angle", default=0.005, type=float_or_none,
help="Specify angular width (in Rad)")
+ parser.add_argument("--xpos", default=None, type=float,
+ help="Specify x position of particles (in mm)")
+ parser.add_argument("--ypos", default=None, type=float,
+ help="Specify y position of particles (in mm)")
parser.add_argument("--zpos", default=None, type=float,
help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index 169badecf..7ec05e2aa 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -148,9 +148,15 @@ if __name__ == '__main__':
# -1000 is safely upstream of detector (to be checked)
# Note zpos is in mm!
- if args.zpos:
+ if args.zpos != None:
sg_dict["z"] = args.zpos
+ if args.xpos != None:
+ sg_dict["x"] = args.xpos
+
+ if args.ypos != None:
+ sg_dict["y"] = args.ypos
+
# Determine energy sampling
if args.sampler == "lin":
sg_dict["energy"] = PG.UniformSampler(args.minE*GeV, args.maxE*GeV)
@@ -159,11 +165,28 @@ if __name__ == '__main__':
elif args.sampler == "const":
sg_dict["energy"] = PG.ConstSampler(args.maxE*GeV)
elif args.sampler == "hist":
- fname, hname = args.hist_name.split(":")
- sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+ nargs = len(args.hist_name.split(":"))
+ if nargs == 2:
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+ elif nargs == 3:
+ fname, hname, scale = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH1Sampler(fname, hname, scale)
+ else:
+ print(f"Can't parse histogram {args.hist_name}!")
+ sys.exit(1)
+
elif args.sampler == "hist2D":
- fname, hname = args.hist_name.split(":")
- sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+ nargs = len(args.hist_name.split(":"))
+ if nargs == 2:
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+ elif nargs == 4:
+ fname, hname, scalex, scaley = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH2Sampler(fname, hname, scalex, scaley)
+ else:
+ print(f"Can't parse histogram {args.hist_name}!")
+ sys.exit(1)
else:
print(f"Sampler {args.sampler} not known!")
sys.exit(1)
@@ -255,8 +278,12 @@ if __name__ == '__main__':
b = time.time()
log.info("Run G4FaserAlg in " + str(b-a) + " seconds")
-#
-# Success should be 0
-#
- sys.exit(not sc.isSuccess())
+
+# Signal errors
+if sc.isSuccess():
+ log.info("Execution succeeded")
+ sys.exit(0)
+else:
+ log.info("Execution failed, return 1")
+ sys.exit(1)
diff --git a/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
index 9b8676752..441ef2923 100755
--- a/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
+++ b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
@@ -185,9 +185,12 @@ cd "${config_file_stem}-${seg_str}"
# Run job
if [[ -z "$tag" ]]; then
faser_particlegun.py "--conf=$config_path" "--segment=$seg_str"
+ gen_code=$?
else
faser_particlegun.py "--conf=$config_path" "--segment=$seg_str" "--tag=$tag"
+ gen_code=$?
fi
+echo "Return code: $gen_code"
#
# Print out ending time
date
@@ -198,25 +201,51 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-HITS.root to $outdest"
- mkdir -p $outdest
- eos cp *-HITS.root ${outdest}/ || true
+ thefile=`ls *-HITS.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that it worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-xAOD.root returned nothing!"
+ copy_code=1
+ fi
+
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
- mkdir -p $logdest
+ eos mkdir -p $logdest
eos cp $logfile $logdest/$logfile
elif ! [ -z "$outdest" ]
then
cd ..
ls -l
echo "copy $logfile to $outdest"
- mkdir -p $outdest
+ eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $gen_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
index 9f541da7c..10549ef0b 100644
--- a/Generators/ParticleGun/python/histsampling.py
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -14,6 +14,7 @@ def load_hist(*args):
"""
Load a histogram from a filename/TFile and histo name. If a single arg is
provided, it has to be a histo object and will be cloned before return.
+
"""
h = None
if len(args) == 1 and issubclass(type(args[0]), ROOT.TH1):
@@ -25,6 +26,7 @@ def load_hist(*args):
#f.Close()
elif type(args[0]) is ROOT.TFile and type(args[1]) is str:
h = args[0].Get(args[1]).Clone()
+
if h is None:
raise Exception("Error in histogram loading from " + args)
else:
@@ -77,7 +79,7 @@ def get_random_bin(globalbins, cheights):
raise Exception("Sample fell outside range of cumulative distribution?!?!")
-def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
+def get_random_x(h, globalbins, cheights, globalbin_to_axisbin, scale=1.):
"""
Choose a random bin via get_random_bin, then pick a uniform random x
point in that bin (without any attempt at estimating the in-bin distribution).
@@ -86,10 +88,10 @@ def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
axisids = globalbin_to_axisbin.get(irand)
assert axisids is not None
xrand = random.uniform(h.GetXaxis().GetBinLowEdge(axisids[0]), h.GetXaxis().GetBinUpEdge(axisids[0]))
- return xrand
+ return scale * xrand
-def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
+def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin, xscale=1., yscale=1.):
"""
Choose a random bin via get_random_bin, then pick a uniform random x,y
point in that bin (without any attempt at estimating the in-bin distribution).
@@ -99,34 +101,71 @@ def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
assert axisids is not None
xrand = random.uniform(h2.GetXaxis().GetBinLowEdge(axisids[0]), h2.GetXaxis().GetBinUpEdge(axisids[0]))
yrand = random.uniform(h2.GetYaxis().GetBinLowEdge(axisids[1]), h2.GetYaxis().GetBinUpEdge(axisids[1]))
- return xrand, yrand
+ return (xscale*xrand), (yscale*yrand)
class TH1(object):
"Minimal wrapper for ROOT TH1, for sampling consistency and easy loading"
def __init__(self, *args):
- self.th1 = load_hist(*args)
+ """ Args can be variable length, but these are now kwargs, so order matters
+
+ fname - filename
+ hname - histogram name
+ xscale - scaling factor for x axis variable
+ """
+ if len(args) > 2:
+ self.th1 = load_hist(*(args[0:2]))
+ else:
+ self.th1 = load_hist(*args)
+
self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+ if len(args) >= 3:
+ self.xscale = float(args[2])
+ else:
+ self.xscale = 1.
+
def GetRandom(self):
"A GetRandom that works for TH1s and uses Python random numbers"
if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
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)
+ return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale)
class TH2(object):
"Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
def __init__(self, *args):
- self.th2 = load_hist(*args)
+ """ Args can be variable length, but these are now kwargs, so order matters
+
+ fname - filename
+ hname - histogram name
+ xscale - scaling factor for x axis variable
+ yscale - scaling factor for y axis variable
+ """
+ if len(args) > 2:
+ self.th2 = load_hist(*(args[0:2]))
+ else:
+ self.th2 = load_hist(*args)
+
self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+ if len(args) >= 3:
+ self.xscale = float(args[2])
+ else:
+ self.xscale = 1.
+
+ if len(args) >= 4:
+ self.yscale = float(args[3])
+ else:
+ self.yscale = 1.
+
+
def GetRandom(self):
"A GetRandom that works for TH2s"
if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
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)
+ return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale, self.yscale)
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index d6e540ea6..9126b607c 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -633,14 +633,23 @@ class EThetaMPhiSampler(MomSampler):
pt = p sin(theta)
"""
- 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 )
+
+ count = 0
+ e = -1
+ while (e < m and count < 5):
+ count += 1
+ if self._theta is None:
+ e,theta = self.energy()
+ else:
+ e = self.energy()
+ theta = self.theta()
+
+ try:
+ p = math.sqrt( e**2 - m**2 )
+ except Exception:
+ raise Exception(f"Error generating E: {e} m: {m}!")
+
pz = p * math.cos(theta)
pt = p * math.sin(theta)
phi = self.phi()
--
GitLab