Update foresee sample config to 13.6 TeV and add ability to sample hists in 1 or 2D

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,

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")

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
 

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", "")

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)
+

Generators/ParticleGun/python/samplers.py

 # 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()