diff --git a/Generators/DIFGenerator/python/DIFSampler.py b/Generators/DIFGenerator/python/DIFSampler.py
index 24659c0178d8cbefc9843a24a1a10ef38bbcde21..69fe82660a4b9d9ae555799d2a8f76dbca5629c6 100644
--- a/Generators/DIFGenerator/python/DIFSampler.py
+++ b/Generators/DIFGenerator/python/DIFSampler.py
@@ -111,7 +111,7 @@ class DIFSampler(PG.ParticleSampler):
return sqrt(m0**4 - (2*m0**2 * m1**2) + (m1**4) - (2*m0**2 * m2**2) - (2*m1**2 *m2**2) + (m2**4)) / (2*m0)
def lorentz_transformation(self):
- self.mother_mom = self.mother.mom.shoot()
+ #self.mother_mom = self.mother.mom.shoot()
Bx = self.mother_mom.Px() / self.mother_mom.E()
By = self.mother_mom.Py() / self.mother_mom.E()
Bz = self.mother_mom.Pz() / self.mother_mom.E()
@@ -132,6 +132,11 @@ class DIFSampler(PG.ParticleSampler):
## The magnitude of the momentum will be equal and opposite
self.mother = self.mother_sampler
+ self.mother.mass_override = False
+ mother_part = self.mother.shoot()[0]
+ self.mother_mom = mother_part.mom
+ self.mother_pos = mother_part.pos
+
p = self.calculate_decay_p(self.mother.mom.mass(),self.daughter1.mass,self.daughter2.mass)
self.daughter1.E = sqrt(self.daughter1.mass**2 + p**2)
@@ -412,4 +417,4 @@ if __name__ == "__main__":
DIFS = DIFSampler()
DIFS.mother_sampler = SP()
- print("DIFSampler: All unit tests passed")
\ No newline at end of file
+ print("DIFSampler: All unit tests passed")
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index 086d3c6aa5e3b5827df893a727d7e6a59f858780..6357163705cd12872a421651495c32db58f021c7 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -159,30 +159,25 @@ def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
pg = cfg.getPrimary()
- from DIFGenerator.DIFSampler import CylinderSampler
-
- from DIFGenerator import DIFSampler
- pg.sampler = DIFSampler(
- daughter1_pid = kwargs.get("daughter1_pid", 11),
- daughter2_pid = kwargs.get("daughter2_pid", -11),
- mother_pid = kwargs.get("mother_pid", None),
- mother_pos = kwargs.get("mother_pos", CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0))
- )
-
from ForeseeGenerator.ForeseeSampler import ForeseeNumpySampler
- mother_mom = ForeseeNumpySampler(
+ mother_part = ForeseeNumpySampler(
model_path = kwargs.get("model_path", "."),
model_name = kwargs.get("model_name", "DarkPhoton"),
com_energy = kwargs.get("com_energy", "14"),
mother_mass = kwargs.get("mother_mass", 0.01),
+ mother_pid = kwargs.get("mother_pid", None),
daughter1_pid = kwargs.get("daughter1_pid", 11),
daughter2_pid = kwargs.get("daughter2_pid", -11),
randomSeed = kwargs.get("randomSeed", None)
)
-
- pg.sampler.mother_sampler = PG.ParticleSampler(pid=kwargs.get("mother_pid", None),mom=mother_mom,n=1,pos=pg.sampler.mother_sampler.pos)
- pg.sampler.mother_sampler.mass_override = False
+ from DIFGenerator import DIFSampler
+ pg.sampler = DIFSampler(
+ daughter1_pid = kwargs.get("daughter1_pid", 11),
+ daughter2_pid = kwargs.get("daughter2_pid", -11),
+ )
+
+ pg.sampler.mother_sampler = mother_part
return cfg
diff --git a/Generators/FlukaReader/CMakeLists.txt b/Generators/FlukaReader/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eebebd79e2b466203bb35ba1aa3d3d265d39920d
--- /dev/null
+++ b/Generators/FlukaReader/CMakeLists.txt
@@ -0,0 +1,11 @@
+################################################################################
+# Package: FlukaReader
+################################################################################
+
+# Declare the package name:
+atlas_subdir( FlukaReader )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+#atlas_install_joboptions( share/*.py )
diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py
new file mode 100644
index 0000000000000000000000000000000000000000..429a443d9bacad678356d794054018a59ea815aa
--- /dev/null
+++ b/Generators/FlukaReader/python/FlukaReaderAlg.py
@@ -0,0 +1,395 @@
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+from GeneratorModules.EvgenAlg import EvgenAlg
+from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType
+from AthenaCommon.SystemOfUnits import GeV, MeV, cm
+from AthenaCommon.Constants import DEBUG
+
+from FaserCosmicGenerator import Range
+
+import ROOT
+
+import numpy as np
+import math
+
+class FlukaReader(EvgenAlg):
+ def __init__(self, name="FlukaReader", MCEventKey="BeamTruthEvent", file_name = "", dist = 0, randomSeed = None, nsamples = 1, test = False):
+ super(FlukaReader,self).__init__(name=name)
+ self.McEventKey = MCEventKey
+ self.file_name = file_name
+ self.dist = dist * 100
+ self.isample = 0
+ self.nsamples = nsamples
+ self.test = test
+
+ self.columns = ["run", "event", "type", "gen", "E", "w", "x", "y", "cosX", "cosY", "age", "_"]
+
+ if randomSeed is not None:
+ self.msg.info(f"Setting seed to {randomSeed}")
+ self.rng = np.random.default_rng(randomSeed)
+ else:
+ self.rng = np.random.default_rng()
+
+ self.file = open(self.file_name)
+
+ if self.test:
+ self.before = dict(zip(self.columns, [[] for i in range(len(self.columns))]))
+ self.after = dict(zip(self.columns, [[] for i in range(len(self.columns))]))
+
+ return
+
+ def genFinalize(self):
+
+ if self.test:
+ self.plot()
+
+ self.file.close()
+ return StatusCode.Success
+
+ def fillEvent(self, evt):
+ "This is called for every real event * the number of samplings"
+
+ # If the sample gets to the number requested, then reset to 0
+ if self.isample == self.nsamples:
+ self.msg.debug("Reseting samples")
+ self.isample = 0
+
+ # Only if the sample is 0 load the new fluka entry
+ if self.isample == 0:
+ self.msg.debug("Loading new fluka event")
+ try:
+ l = next(self.file)
+ except StopIteration:
+ return StatusCode.Success
+
+ entry = dict(zip(self.columns, l.strip("\n").split()))
+ for i,c in enumerate(self.columns):
+ if i < 4:
+ entry[c] = int(entry[c])
+ else:
+ entry[c] = float(entry[c])
+
+ self.entry = entry
+
+ # Call for each sample of each event
+ self.msg.debug(f"Processing sample {self.isample}")
+ self.process(self.entry, evt)
+ self.isample += 1
+
+ return StatusCode.Success
+
+ def angle(self, cosTheta):
+ "Convert cos(theta) wrt x or y axis to theta wrt z axis"
+ return np.pi/2 - np.arccos(cosTheta)
+
+ def pid(self, ftype):
+ "Convert fluka particle type to PID"
+ if ftype == 10: # mu+
+ return -13
+ elif ftype == 11: # mu-
+ return 13
+ else:
+ return 0
+
+ def path_length(self, z, cosThetaX, cosThetaY):
+ "Get path length traversed in the material, taking into account incident angle"
+
+ # Convert theta wrt x and y axis to wrt z axis
+ thetaX = self.angle(cosThetaX)
+ thetaY = self.angle(cosThetaY)
+
+ # Correct z for angle
+ zcorr = z / np.abs(np.cos(thetaX)) / np.abs(np.cos(thetaY))
+
+ return zcorr
+
+ def energy_after_loss_exact(self, e, zcorr):
+ "Calculate exact energy after loss in material"
+ return Range.muPropagate(e, zcorr/1000.)
+
+ def energy_after_loss(self, e, cosThetaX, cosThetaY, zcorr, a = 2e-3, b = 4e-6):
+ "Calculate approximate energy after loss in material"
+
+ # Based on
+ # http://www.bartol.udel.edu/~stanev/lectures/apr17.pdf
+ # and PDG chapter 27 fig 27.1
+ # https://www.google.co.uk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiG9YjNtvr2AhUbiVwKHdQfD9sQFnoECAsQAQ&url=https%3A%2F%2Fpdg.lbl.gov%2F2005%2Freviews%2Fpassagerpp.pdf&usg=AOvVaw1HGA5PZtC2UiqA6B7_C5dz
+
+ eps = a/b
+ return (e + eps) * np.exp(-b * zcorr) - eps
+
+ def mean_scattering_angle(self, e, cosThetaX, cosThetaY, zcorr, X0 = 10.02, m = 105.66e-3, charge = 1, beta = 1):
+ "Calculate mean scattering angle over many scatters for given energy and length z"
+
+ # Based on PDG chapter 27 eqns 27.10, 27.16, 27.17
+ # https://www.google.co.uk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiG9YjNtvr2AhUbiVwKHdQfD9sQFnoECAsQAQ&url=https%3A%2F%2Fpdg.lbl.gov%2F2005%2Freviews%2Fpassagerpp.pdf&usg=AOvVaw1HGA5PZtC2UiqA6B7_C5dz
+ # and
+ # https://pdg.lbl.gov/2014/AtomicNuclearProperties/HTML/standard_rock.html
+
+ # Convert E to momentum [GeV]
+ p = np.sqrt(e**2 - m**2)
+
+ # Mean angle [GeV and cm]
+ c = 1 # n.u
+ theta0 = 13.6e-3 / (p * c * beta) * charge * np.sqrt(zcorr/X0) * (1 + 0.38 * math.log(zcorr, X0))
+ return theta0
+
+ def scattering_angle(self, cosTheta, theta0, rand1):
+ "Calculate actual scattering angle over many scatters for given start angle and mean angle"
+
+ # Convert theta wrt x or y axis to wrt z axis
+ theta = self.angle(cosTheta)
+
+ # Add random scattering angle
+ return theta + rand1 * theta0
+
+ def scattering_postition(self, x, cosThetaX, cosThetaY, zcorr, theta0, rand1, rand2):
+ "Calculate transverse scattering position over many scatters for given start angle and mean angle + length z"
+
+ # Convert theta wrt x to wrt z axis
+ thetaX = self.angle(cosThetaX)
+
+ xout = np.copy(x)
+
+ # Add displacement due to initial angle
+ #xang = z * np.tan(thetaX)
+ xang = zcorr * np.sin(thetaX)
+ xout += xang
+
+ # Add displacement due to multiple scattering
+ dx = rand1 * zcorr * theta0 / np.sqrt(12) + rand2 * zcorr * theta0/2
+ xout += dx
+
+ return xout
+
+ def propagate(self, entry):
+ "Propagate the particle through a given distance of standard rock using the small-angle approxiumation"
+
+ if self.dist == 0:
+ return entry
+
+ # Random numbers
+ rand1 = self.rng.normal(0, 1)
+ rand2 = self.rng.normal(0, 1)
+
+ # Get entry info
+ e = entry["E"]
+ x = entry["x"]
+ y = entry["y"]
+ cosX = entry["cosX"]
+ cosY = entry["cosY"]
+
+ # Correct path length for angles
+ z = self.path_length(self.dist, cosX, cosY)
+
+ # Account for energy loss
+ #eout = self.energy_after_loss(e, cosX, cosY, z)
+ eout = self.energy_after_loss_exact(e, z)
+
+ # Account for scattering on angle and position
+ theta0 = self.mean_scattering_angle(e, cosX, cosY, z)
+ thetaXout = self.scattering_angle(cosX, theta0, rand1)
+ thetaYout = self.scattering_angle(cosY, theta0, rand1)
+ xout = self.scattering_postition(x, cosX, cosY, z, theta0, rand1, rand2)
+ yout = self.scattering_postition(y, cosY, cosX, z, theta0, rand1, rand2)
+
+ # Update entry info
+ entry["E"] = eout
+ entry["x"] = xout
+ entry["y"] = yout
+ entry["cosX"] = np.cos(np.pi/2 + thetaXout)
+ entry["cosY"] = np.cos(np.pi/2 + thetaYout)
+
+ return entry
+
+
+ def process(self, entry, evt):
+
+ if self.test:
+ for k,v in entry.items():
+ self.before[k].append(float(v))
+
+ if self.msg.level > DEBUG:
+ print("Original Entry", entry)
+
+ if self.dist != 0:
+ # Propoagate to FASER
+ entry = self.propagate(entry)
+ elif self.nsamples != 1:
+ # Else smear if sampling more than once
+ entry["E"] *= self.rng.normal(1, 0.05)
+ entry["x"] *= self.rng.normal(1, 0.05)
+ entry["y"] *= self.rng.normal(1, 0.05)
+ entry["cosX"] = np.cos(np.arccos(entry["cosX"]) * self.rng.normal(1, 0.05))
+ entry["cosY"] = np.cos(np.arccos(entry["cosY"]) * self.rng.normal(1, 0.05))
+
+ if self.msg.level > DEBUG:
+ print("Propagated/Smeared Entry", entry)
+
+
+ if self.test:
+ for k,v in entry.items():
+ self.after[k].append(float(v))
+
+ try:
+ from AthenaPython.PyAthena import HepMC3 as HepMC
+ except ImportError:
+ from AthenaPython.PyAthena import HepMC as HepMC
+
+ # Add weight, correcting for mutliple sampling
+ evt.weights().push_back(entry["w"] / self.nsamples)
+
+
+ # Setup MC event
+ mcEventType = EventType()
+ mcEventType.add_type(EventType.IS_SIMULATION)
+
+ mcEventId = EventID(run_number = entry["run"], event_number = entry["event"])
+
+ mcEventInfo = EventInfo(id = mcEventId, type = mcEventType)
+
+ self.evtStore.record(mcEventInfo, "McEventInfo", True, False)
+
+ ROOT.SetOwnership(mcEventType, False)
+ ROOT.SetOwnership(mcEventId, False)
+ ROOT.SetOwnership(mcEventInfo, False)
+
+ # Create HepMC Vertex
+ z = 0
+ pos = HepMC.FourVector(entry["x"] * cm, entry["y"] * cm, z, 0)
+ gv = HepMC.GenVertex(pos)
+
+ ROOT.SetOwnership(gv, False)
+ evt.add_vertex(gv)
+
+ # Create HepMC particle
+ gp = HepMC.GenParticle()
+
+ m = 105.66
+ e = entry["E"] * 1000.
+ p = np.sqrt(e**2 - m**2)
+
+ thetaX = self.angle(entry["cosX"])
+ thetaY = self.angle(entry["cosY"])
+ px = p * np.sin(thetaX) * np.cos(thetaY)
+ py = p * np.cos(thetaX) * np.sin(thetaY)
+ pz = p * np.cos(thetaX) * np.cos(thetaY)
+
+ mom = HepMC.FourVector(px, py, pz, e)
+
+ gp.set_momentum(mom)
+ gp.set_generated_mass(m)
+ gp.set_pdg_id(self.pid(entry["type"]))
+ gp.set_status(1)
+
+ ROOT.SetOwnership(gp, False)
+ gv.add_particle_out(gp)
+
+ return
+
+ def plot(self):
+ "Plot entries before and after propagation/smeating for tests"
+
+ if not self.test:
+ return
+
+ import matplotlib.pyplot as plt
+
+ plt.figure()
+ ebins = np.linspace(0, 5000, 50)
+ plt.hist(self.before["E"], bins=ebins, histtype='step', color = "g", fill = False)
+ plt.hist(self.after["E"], bins = ebins, histtype='step', color = "r", fill = False)
+ plt.gca().set_yscale('log')
+ plt.savefig("energy.eps")
+
+ plt.figure()
+ 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)
+ plt.hist(thetaXout, bins = tbins, histtype='step', color = "r", fill = False)
+ plt.gca().set_yscale('log')
+ plt.savefig("thetaX.eps")
+
+ plt.figure()
+ 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)
+ plt.hist(thetaYout, bins = tbins, histtype='step', color = "r", fill = False)
+ plt.gca().set_yscale('log')
+ plt.savefig("thetaY.eps")
+
+ plt.figure()
+ xbins = np.linspace(-300, 300, 100)
+ plt.hist(self.before["x"], bins = xbins, histtype='step', color = "g", fill = False)
+ plt.hist(self.after["x"], bins = xbins, histtype='step', color = "r", fill = False)
+ plt.gca().set_yscale('log')
+ plt.savefig("x.eps")
+
+ plt.figure()
+ plt.hist(self.before["y"], bins = xbins, histtype='step', color = "g", fill = False)
+ plt.hist(self.after["y"], bins = xbins, histtype='step', color = "r", fill = False)
+ plt.gca().set_yscale('log')
+ plt.savefig("y.eps")
+
+ return
+
+if __name__ == "__main__":
+
+# from AthenaCommon.AlgSequence import AlgSequence
+# job = AlgSequence()
+# job += FlukaReader(file_name = "/user/gwilliam/unit30_Nm", dist = (480-409))
+#
+# from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
+# ostream = AthenaPoolOutputStream( "StreamEVGEN" , "evgen.pool.root", noTag=True )
+# ostream.ItemList.remove("EventInfo#*")
+# ostream.ItemList += [ "EventInfo#McEventInfo",
+# "McEventCollection#*" ]
+#
+# theApp.EvtMax = 1000
+
+
+ import argparse, sys
+ parser = argparse.ArgumentParser(description="Run Fluka reader")
+ parser.add_argument("file", help = "Path to fluka file")
+ parser.add_argument("--dist", "-z", default = 0, type = float, help = "depth of standard rock to propagate through [m]")
+ parser.add_argument("--output", "-o", default = "evgen.pool.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")
+ parser.add_argument("--randomSeed", "-r", default=12345, type=int, help = "Seed for random number generator")
+ parser.add_argument("--nsamples", "-s", default = 1, type = int, help = "Number of times to sample each event")
+ parser.add_argument("--test", "-t", action = "store_true", help = "Make test plots")
+ args = parser.parse_args()
+
+
+ from AthenaCommon.Logging import log
+ from AthenaCommon.Constants import DEBUG, INFO
+
+ from AthenaCommon.Configurable import Configurable
+ Configurable.configurableRun3Behavior = 1
+
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ 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 = args.output
+ ConfigFlags.lock()
+
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ cfg = MainServicesCfg(ConfigFlags)
+
+ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+ from AthenaConfiguration.ComponentFactory import CompFactory
+
+ acc = ComponentAccumulator()
+ reader = FlukaReader("FlukaReader", MCEventKey=args.mcEventKey, file_name = args.file, dist = args.dist, randomSeed = args.randomSeed, nsamples = args.nsamples, test = args.test)
+ reader.OutputLevel = INFO
+ acc.addEventAlgo(reader)
+ cfg.merge(acc)
+
+ itemList = [ "EventInfo#McEventInfo", "McEventCollection#*" ]
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ cfg.merge(OutputStreamCfg(ConfigFlags, "EVNT", itemList, disableEventTag = True))
+
+ sc = cfg.run(maxEvents = args.nevents * args.nsamples)
+ sys.exit(not sc.isSuccess())
diff --git a/Generators/ForeseeGenerator/python/ForeseeSampler.py b/Generators/ForeseeGenerator/python/ForeseeSampler.py
index d80be0aa3aa080d93eb9ec7510b61aebf6089464..2f7d7241e560e552b8df3bcd9cb9b62e321f17ff 100644
--- a/Generators/ForeseeGenerator/python/ForeseeSampler.py
+++ b/Generators/ForeseeGenerator/python/ForeseeSampler.py
@@ -4,18 +4,24 @@ import numpy as np
import ParticleGun as PG
-class ForeseeNumpySampler(PG.MomSampler):
+from DIFGenerator.DIFSampler import CylinderSampler
+
+class ForeseeNumpySampler(PG.ParticleSampler):
""" Sample from the output of Foresee generation in numpy format with columns E, theta, weight"""
def __init__(self, model_path = ".", model_name = "DarkPhoton", com_energy = "14", mother_mass = 0.01,
- coupling = None, daughter1_pid = 11, daughter2_pid = -11, randomSeed = None):
+ coupling = None, mother_pid = None, daughter1_pid = 11, daughter2_pid = -11, randomSeed = None):
self.path = os.path.expanduser(os.path.expandvars(model_path))
self.modelname = model_name
self.energy = com_energy
- self._mass = mother_mass
+ self.mass = mother_mass
self.coupling = coupling
+ self.pid = mother_pid
self.daughter1_pid = daughter1_pid
self.daughter2_pid = daughter2_pid
+ self.n = 1
+ self.distance = 480 # m
+ self.mass_override = False
if randomSeed is not None:
print(f"Setting seed to {randomSeed}")
@@ -33,9 +39,9 @@ class ForeseeNumpySampler(PG.MomSampler):
if self.path.endswith(".npy"):
filename = self.path
elif self.coupling is None:
- filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self._mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+ filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self.mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
else:
- filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self._mass}GeV_c{self.coupling}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+ filename = f"{self.path}/files/models/{self.modelname}/events/events_{self.energy}TeV_m{self.mass}GeV_c{self.coupling}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
print(f"Reading data from file: {filename}")
self.data = np.load(filename)
@@ -45,22 +51,42 @@ class ForeseeNumpySampler(PG.MomSampler):
return
- def mass(self):
- "Mass converted from GeV to MeV"
- return self._mass * 1000
+# def mass(self):
+# "Mass converted from GeV to MeV"
+# return self._mass * 1000
def shoot(self):
"Choose a random item from the data, based on the probability"
- # TODO: what about reuse of events?
energy, theta, weight = self.rng.choice(self.data, axis = 1, p = self.prob)
self.xs += weight
+ # Convert mass to MeV
+ mass = self.mass * 1000
+
+ # Sample phi
+ phi = self.rng.uniform(0, 2*np.pi)
+
+ # Sample z
+ z = self.rng.uniform(-1500, 0)
+
# Create mother momentum, randomly sampling phi
- mother_mom = PG.EThetaMPhiSampler(energy, theta, self.mass(), [0, 2*np.pi])
+ self.mom = PG.EThetaMPhiSampler(energy, theta, mass, phi)
- return mother_mom.shoot()
+ # Create mother pos, randomly sampling phi
+ r = (self.distance * 1000 + abs(z)) * np.tan(theta)
+
+ self.pos = CylinderSampler(r**2, phi, z, 0)
+
+ # Create particle
+ p = PG.SampledParticle(self.pid)
+ p.mom = self.mom.shoot()
+ p.pos = self.pos.shoot()
+ p.mass = mass
+ #self.mom.mass = mass
+
+ return [p]
class ForeseeSampler(PG.MomSampler):
"""Create events from foresee directly on the fly
@@ -266,7 +292,7 @@ if __name__ == "__main__":
path = os.path.expandvars("$Calypso_DIR/../calypso/Generators/ForeseeGenerator/data/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy")
-
+ path = "files/models/DarkPhoton/events/events_14TeV_m0.1GeV_c1e-05to_11_-11.npy"
modelname = "DarkPhoton"
mass = 0.1
@@ -281,21 +307,17 @@ if __name__ == "__main__":
d1mom = []
# Accounting for rounding
- epsilon = 5
+ epsilon = 6
- # Mother position sampler over Z
- mother_pos = PG.PosSampler(x=0,y=0,z=[-1500,0],t=0)
-
- # Create mother momentum sampler reading data from foresee
- mother_mom_sampler = ForeseeNumpySampler(model_path = path, model_name = modelname, com_energy = "14", mother_mass = 0.1, coupling = 1e-5, daughter1_pid = 11, daughter2_pid = -11)
+ # Create mother sampler reading data from foresee
+ mother_sampler = ForeseeNumpySampler(model_path = path, model_name = modelname, com_energy = "14", mother_mass = 0.1, coupling = 1e-5, mother_pid = None, daughter1_pid = 11, daughter2_pid = -11)
# Create decay-in-flight
d = DIFSampler(11, -11, None)
- d.mother_sampler = PG.ParticleSampler(pid=None,mom=mother_mom_sampler,n=1,pos=mother_pos)
- d.mother_sampler.mass_override = False
+ d.mother_sampler = mother_sampler
# Loop over a range of events
- for i in range(1000):
+ for i in range(100000):
# Shoot the decay in flight
daughters = d.shoot()
@@ -434,7 +456,7 @@ if __name__ == "__main__":
plt.savefig(f"{modelname}_PG_d1_m{mass}.png")
- print (f"x-sect = {mother_mom_sampler.xs} pb")
+ print (f"x-sect = {mother_sampler.xs} pb")
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 8b53c84caad9190e4a9bda75f05970ec5f524da8..9c6be63475b6e066ece1ceac0aea059b8837a0dc 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -63,11 +63,11 @@ class HistSvc(object):
class EvgenValidation(EvgenAnalysisAlg):
"Gen-level validation"
- def __init__(self, name = "EvgenValidation"):
+ def __init__(self, name = "EvgenValidation", ndaughters = 2, outname = "validation.root"):
super(EvgenValidation, self).__init__(name=name)
self.hists = HistSvc()
- self.ndaughter = 2
- self.outname = "validation.root"
+ self.ndaughters = ndaughters
+ self.outname = outname
def binning(self):
"binning for theta vs phi plot"
@@ -84,8 +84,10 @@ class EvgenValidation(EvgenAnalysisAlg):
# Daughter i
tbins, pbins = self.binning()
- for i in range(self.ndaughter):
+ for i in range(self.ndaughters):
+ self.hists.add(f"E_d{i}", 100, 0, 10000)
self.hists.add(f"P_d{i}", 100, 0, 10000)
+ self.hists.add(f"Pz_d{i}", 100, 0, 10000)
self.hists.add(f"Pt_d{i}", 100, 0, 1)
self.hists.add(f"Theta_d{i}", 20, 0, 0.001)
self.hists.add(f"Phi_d{i}", 16, -3.2, 3.2)
@@ -93,53 +95,82 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists.add(f"Mass_d{i}", 200, 0, 0.01)
# Mother
- self.hists.add("P_M", 100, 0, 10000)
+ 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("Phi_M", 16, -3.2, 3.2)
self.hists.add("Mass_M", 200, 0, 1)
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)
return StatusCode.Success
def fillKin(self, label, p, mass = True, twoD = True):
- self.hists[f"P_{label}"].Fill(p.rho()/1000, 1)
- self.hists[f"Pt_{label}"].Fill(p.perp()/1000, 1)
- self.hists[f"Theta_{label}"].Fill(p.theta(), 1)
- self.hists[f"Phi_{label}"].Fill(p.phi(), 1)
+ self.hists[f"E_{label}"].Fill(p.t()/1000, self.weight)
+ self.hists[f"P_{label}"].Fill(p.rho()/1000, self.weight)
+ self.hists[f"Pz_{label}"].Fill(p.pz()/1000, self.weight)
+ self.hists[f"Pt_{label}"].Fill(p.perp()/1000, self.weight)
+ self.hists[f"Theta_{label}"].Fill(p.theta(), self.weight)
+ self.hists[f"Phi_{label}"].Fill(p.phi(), self.weight)
if mass:
- self.hists[f"Mass_{label}"].Fill(p.m()/1000, 1)
+ self.hists[f"Mass_{label}"].Fill(p.m()/1000, self.weight)
if twoD:
- self.hists[f"ThetaVsP_{label}"].Fill(p.theta(), p.rho()/1000, 1)
+ self.hists[f"ThetaVsP_{label}"].Fill(p.theta(), p.rho()/1000, self.weight)
return
def fillDaughter(self, p):
- self.hists["PIDs"].Fill(p.pdg_id())
+ 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)
+ return
+
+
def execute(self):
evt = self.events()[0]
+ self.weight = evt.weights()[0]
# Loop over all particles in events (assuming mother not stored)
momenta = []
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:
+ #p.production_vertex().print()
+ llp_vtx = p.production_vertex().point3d()
# Fill daughter plots
- for i in range(self.ndaughter):
+ for i in range(self.ndaughters):
self.fillKin(f"d{i}", momenta[i])
# Fill mother plots
self.fillKin("M", mother, mass = True)
+
+ # Fill vertex plots
+ self.fillVertex(llp_vtx)
return StatusCode.Success
@@ -153,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 = float, help = "Number of daugthers to plot")
+ parser.add_argument("--ndaugthers", "-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")
@@ -187,10 +218,8 @@ if __name__ == "__main__":
fix()
acc = ComponentAccumulator()
- valid = EvgenValidation()
+ valid = EvgenValidation("EvgenValidation", ndaughters = args.ndaugthers, outname = args.output)
valid.McEventKey = args.mcEventKey
- valid.ndaughters = args.ndaugthers
- valid.outname = args.output
acc.addEventAlgo(valid)
cfg.merge(acc)
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 2a9e7a7fefd775c9d012fdfefe38c4ca8227b6d7..ba63087bf08bf1e2571865651f255cd9ad177084 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -39,7 +39,9 @@ class ForeseeGenerator(object):
# Set detector ...
self.foresee = Foresee()
- self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1", channels=[self.mode], distance=480, length=1.5 , luminosity=150)
+ self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1",
+ channels=[self.mode], distance=480, length=1.5 ,
+ luminosity=1/1000.) # 1 pb-1
# Set model ...
self.model = Model(self.modelname)
@@ -72,7 +74,6 @@ class ForeseeGenerator(object):
energy = self.energy,
nsample = 10)
- # This is a bit handwavey according to Felix
self.model.add_production_mixing(
pid = "113",
mixing = "coupling * 0.3/5. * 0.77545**2/abs(mass**2-0.77545**2+0.77545*0.147*1j)",
@@ -80,7 +81,6 @@ class ForeseeGenerator(object):
energy = self.energy,
)
- # Question on validity as FASER gets larger but OK for currrent size according to Felix
self.model.add_production_direct(
label = "Brem",
energy = self.energy,
@@ -153,15 +153,17 @@ class ForeseeGenerator(object):
)
# Get LLP spectrum
- self.foresee.set_model(model=self.model)
- plt = self.foresee.get_llp_spectrum(self.mass, coupling=1, do_plot=True) # This is just a reference coupling
+ self.foresee.set_model(model=self.model)
+ # This is just a reference coupling
+ plt = self.foresee.get_llp_spectrum(self.mass, coupling=1, do_plot=True)
plt.savefig(f"{self.modelname}_m{self.mass}.png")
def flatten(l):
return [i for sublist in l for i in sublist]
# Get list of events within detector
- coups, ctaus, nsigs, energies, weights, thetas = self.foresee.get_events(mass=self.mass, energy=self.energy, couplings=self.couplings)
+ output = self.foresee.get_events(mass=self.mass, energy=self.energy, couplings=self.couplings)
+ coups, ctaus, nsigs, energies, weights, thetas = output
self.plot(flatten(thetas), flatten(energies), flatten(weights))
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index 5b3d8b56f2ec4dc60f308be26374191ca47bfc2e..8d5a04701c82deedadada27fc76d11fa4a6a8d30 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -14,8 +14,10 @@ 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")
-
- if r != 1:
+
+ if isinstance(r, tuple):
+ h.Rebin2D(r[0], r[1])
+ elif r != 1:
h.Rebin(r)
if xtitle is not None:
@@ -110,3 +112,11 @@ if __name__ == "__main__":
]
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),
+ Hist("Vtx_XY", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5))
+ ]
+
+ plotn(f, config, 2, 2, "vtx")
diff --git a/Generators/ParticleGun/CMakeLists.txt b/Generators/ParticleGun/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1794305ed0291e77478aae72b1216db87fce0810
--- /dev/null
+++ b/Generators/ParticleGun/CMakeLists.txt
@@ -0,0 +1,14 @@
+################################################################################
+# Package: ParticleGun
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ParticleGun )
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+
+# Install files from the package:
+atlas_install_joboptions( share/common/*.py
+ share/example/*.py )
+
diff --git a/Generators/ParticleGun/README b/Generators/ParticleGun/README
new file mode 100644
index 0000000000000000000000000000000000000000..6cdb698eda4f549c076ac3a3ce732ad65225c11d
--- /dev/null
+++ b/Generators/ParticleGun/README
@@ -0,0 +1,5 @@
+ParticleGun documentation
+-------------------------
+See https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/ParticleGunForAtlas
+for some coherent documentation that should be kept up to date.
+
diff --git a/Generators/ParticleGun/python/__init__.py b/Generators/ParticleGun/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..4e88fc56c49679a210ad57b62055a6a7c0867ddd
--- /dev/null
+++ b/Generators/ParticleGun/python/__init__.py
@@ -0,0 +1,129 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+from GeneratorModules.EvgenAlg import EvgenAlg
+from ParticleGun.samplers import ParticleSampler
+from ParticleGun.samplers import * # noqa: F401, F403 (import into our namespace)
+# commenting out the HepMC import for now
+#try:
+# from AthenaPython.PyAthena import HepMC3 as HepMC
+#except ImportError:
+# from AthenaPython.PyAthena import HepMC as HepMC
+
+from AthenaPython.PyAthena import StatusCode
+import ROOT,random
+
+__author__ = "Andy Buckley <andy.buckley@cern.ch>"
+
+class ParticleGun(EvgenAlg):
+ """
+ A simple but flexible algorithm for generating events from simple distributions.
+ """
+
+ def __init__(self, name="ParticleGun", randomSvcName="AtRndmGenSvc", randomStream="ParticleGun", randomSeed=None):
+ super(ParticleGun, self).__init__(name=name)
+ self.samplers = [ParticleSampler()]
+ self.randomStream = randomStream
+ self.randomSvcName = randomSvcName
+ self.randomSeed = randomSeed
+
+ @property
+ def sampler(self):
+ "Get the first (and presumed only) sampler"
+ return self.samplers[0] if self.samplers else None
+ @sampler.setter
+ def sampler(self, s):
+ "Set the samplers list to include only a single sampler, s"
+ self.samplers = [s]
+
+
+ def initialize(self):
+ """
+ Pass the AtRndmGenSvc seed to Python's random module, or use a fixed value set via pg.randomSeed.
+ """
+ seed = None
+ ## Use self.randomSeed directly, or if it's None find a seed string from the ATLAS random number service
+ if self.randomSeed is not None:
+ seed = self.randomSeed
+ else:
+ randomSvc = getattr(svcMgr, self.randomSvcName, None)
+ if randomSvc is not None:
+ for seedstr in randomSvc.Seeds:
+ if seedstr.startswith(self.randomStream):
+ seed = seedstr
+ self.msg.info("ParticleGun: Using random seed '%s' ", seed)
+ break
+ if seed is None:
+ self.msg.warning("ParticleGun: Failed to find a seed for the random stream named '%s' ", self.randomStream)
+ else:
+ self.msg.warning("ParticleGun: Failed to find random number service called '%s' ", self.randomSvcName)
+ ## Apply the seed
+ if seed is not None:
+ random.seed(seed)
+ return StatusCode.Success
+ else:
+ self.msg.error("ParticleGun: randomSeed property not set, and no %s random number service found", self.randomSvcName)
+ return StatusCode.Failure
+
+
+ def fillEvent(self, evt):
+ """
+ Sample a list of particle properties, which are then used to create a new GenEvent in StoreGate.
+ """
+ ## Set event weight(s)
+ # TODO: allow weighted sampling?
+ try:
+ from AthenaPython.PyAthena import HepMC3 as HepMC
+ except ImportError:
+ from AthenaPython.PyAthena import HepMC as HepMC
+ evt.weights().push_back(1.0)
+
+ ## Make and fill particles
+ for s in self.samplers:
+ particles = s.shoot()
+ for p in particles:
+ ## Debug printout of particle properties
+ #print DEBUG0, p.pid, p.mom.E(), p.mom.Pt(), p.mom.M()
+ #print "DEBUG1 (px,py,pz,E) = (%0.2e, %0.2e, %0.2e, %0.2e)" % (p.mom.Px(), p.mom.Py(), p.mom.Pz(), p.mom.E())
+ #print "DEBUG2 (eta,phi,pt,m) = (%0.2e, %0.2e, %0.2e, %0.2e)" % (p.mom.Eta(), p.mom.Phi(), p.mom.Pt(), p.mom.M())
+ #print "DEBUG3 (x,y,z,t) = (%0.2e, %0.2e, %0.2e, %0.2e)" % (p.pos.X(), p.pos.Y(), p.pos.Z(), p.pos.T())
+
+ ## Make particle-creation vertex
+ # TODO: do something cleverer than one vertex per particle?
+ pos = HepMC.FourVector(p.pos.X(), p.pos.Y(), p.pos.Z(), p.pos.T())
+ gv = HepMC.GenVertex(pos)
+ ROOT.SetOwnership(gv, False)
+ evt.add_vertex(gv)
+
+ ## Make particle with status == 1
+ mom = HepMC.FourVector(p.mom.Px(), p.mom.Py(), p.mom.Pz(), p.mom.E())
+ gp = HepMC.GenParticle()
+ gp.set_status(1)
+ gp.set_pdg_id(p.pid)
+ gp.set_momentum(mom)
+ if p.mass is not None:
+ gp.set_generated_mass(p.mass)
+ ROOT.SetOwnership(gp, False)
+ gv.add_particle_out(gp)
+
+ return StatusCode.Success
+
+
+## PyAthena HepMC notes
+#
+## evt.print() isn't valid syntax in Python2 due to reserved word
+# TODO: Add a Pythonisation, e.g. evt.py_print()?
+#getattr(evt, 'print')()
+#
+## How to check that the StoreGate key exists and is an McEventCollection
+# if self.sg.contains(McEventCollection, self.sgkey):
+# print self.sgkey + " found!"
+#
+## Modifying an event other than that supplied as an arg
+# mcevts = self.sg[self.sgkey]
+# for vtx in mcevts[0].vertices: # only way to get the first vtx?!
+# gp2 = HepMC.GenParticle()
+# gp2.set_momentum(HepMC.FourVector(1,2,3,4))
+# gp2.set_status(1)
+# vtx.add_particle_out(gp2)
+# break
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
new file mode 100644
index 0000000000000000000000000000000000000000..c64112cc84e23963b63ea381817d1457f1c2a122
--- /dev/null
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -0,0 +1,132 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+"""
+Tools for histogram sampling, in particular inverse transform sampling which is
+missing from ROOT's TH2 classes.
+"""
+
+__author__ = "Andy Buckley <andy.buckley@cern.ch>"
+
+import random, ROOT
+
+
+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):
+ h = args[0].Clone()
+ elif len(args) == 2:
+ if isinstance(args[0], str) and isinstance(args[1], str) :
+ f = ROOT.TFile.Open(args[0])
+ h = f.Get(args[1]).Clone()
+ #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)
+ return h
+
+
+def get_sampling_vars(h):
+ """
+ Get the following from a histogram h, since the ROOT API sucks:
+ * list of global bin IDs (not even contiguous for 2D, gee thanks ROOT)
+ * dict mapping global bin IDs to a tuple of axis bin IDs
+ * list of nbins+1 cumulative bin values, in the same order as globalbins
+ """
+ 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):
+ for ix in range(1, h.GetNbinsX()+1):
+ iglobal = h.GetBin(ix)
+ globalbins.append(iglobal)
+ globalbin_to_axisbin[iglobal] = (ix,)
+ cheights.append(cheights[-1] + h.GetBinContent(iglobal))
+ elif issubclass(type(h), ROOT.TH2):
+ for ix in range(1, h.GetNbinsX()+1):
+ for iy in range(1, h.GetNbinsY()+1):
+ iglobal = h.GetBin(ix, iy)
+ globalbins.append(iglobal)
+ globalbin_to_axisbin[iglobal] = (ix, iy)
+ cheights.append(cheights[-1] + h.GetBinContent(iglobal))
+ return globalbins, globalbin_to_axisbin, cheights
+
+
+def get_random_bin(globalbins, cheights):
+ """
+ Choose a random bin from the cumulative distribution list of nbins+1 entries.
+
+ TODO: Search more efficiently (lin and log guesses, then lin search or
+ binary split depending on vector size).
+ """
+ assert len(cheights) == len(globalbins)+1
+ randomheight = random.uniform(0, cheights[-1])
+ for i, iglobal in enumerate(globalbins):
+ if randomheight >= cheights[i] and randomheight < cheights[i+1]:
+ return iglobal
+ raise Exception("Sample fell outside range of cumulative distribution?!?!")
+
+
+def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
+ """
+ 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).
+ """
+ irand = get_random_bin(globalbins, cheights)
+ 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
+
+
+def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
+ """
+ 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).
+ """
+ irand = get_random_bin(globalbins, cheights)
+ axisids = globalbin_to_axisbin.get(irand)
+ 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
+
+
+class TH1(object):
+ "Minimal wrapper for ROOT TH1, for sampling consistency and easy loading"
+
+ def __init__(self, *args):
+ self.th1 = load_hist(*args)
+ self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+
+ 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)
+
+ 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"
+
+ def __init__(self, *args):
+ self.th2 = load_hist(*args)
+ self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+
+ 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)
+
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..4a232f7d7120136f29f07b493b3dd64c0137d0d3
--- /dev/null
+++ b/Generators/ParticleGun/python/samplers.py
@@ -0,0 +1,909 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+import ROOT, math, random
+from ParticleGun.histsampling import TH1
+
+## For convenience
+PI = math.pi
+TWOPI = 2*math.pi
+
+
+class Sampler(object):
+ "Base class for all samplers"
+
+ def shoot(self):
+ return RuntimeError("Can't sample from an abstract sampler object.")
+
+ def __call__(self):
+ """This is the call method that will actually be used (so that normal
+ functions can also be passed in as samplers)."""
+ return self.shoot()
+
+ # TODO: add a sampling weight?
+
+
+class ConstSampler(Sampler):
+ "A special-case sampler which just returns one value rather than sampling."
+
+ def __init__(self, val):
+ self.val = val
+
+ def shoot(self):
+ return self.val
+
+ def __repr__(self):
+ return "ConstSampler[%s]" % str(self.val)
+
+
+## Continuous distribution samplers
+
+class ContinuousSampler(Sampler):
+ "Base class for samplers from continuous distributions."
+ pass
+
+
+class UniformSampler(ContinuousSampler):
+ "Uniformly sample in the range [low,high)."
+
+ def __init__(self, low, high):
+ assert(low <= high)
+ self.low = float(low)
+ self.high = float(high)
+
+ def shoot(self):
+ return random.uniform(self.low, self.high)
+
+
+class ModUniformSampler(ContinuousSampler):
+ "Uniformly sample in the modulus range (-high,low]+[low,high)."
+
+ def __init__(self, low, high):
+ assert(low == abs(low) and high == abs(high))
+ assert(low <= high)
+ self.low = float(low)
+ self.high = float(high)
+
+ def shoot(self):
+ val = random.uniform(self.low, self.high)
+ if random.random() > 0.5:
+ val *= -1
+ return val
+
+
+class DisjointUniformSampler(ContinuousSampler):
+ "Uniformly sample from a set of disjoint intervals."
+
+ def __init__(self, ranges):
+ """
+ The ranges variable can either be a list of increasing numbers or a
+ list of pairs of numbers.
+
+ The former case will be treated as
+ defining alternating on/off ranges for sampling, starting with an active
+ one (i.e. it's a list of bin edges). The latter way specifically lists
+ the 'on' regions only, with their start and end values in the pairs.
+
+ The behaviour is undefined if the numbers are not ordered or overlap --
+ i.e. it might work but hasn't been designed that way and might change in
+ future. Don't rely on this behaviour!
+ """
+ if not ranges:
+ raise Exception("You must supply at least one non-null sampling range")
+ if hasattr(ranges[0], "__len__"):
+ assert all(len(x) == 2 for x in ranges)
+ self.ranges = ranges
+ else:
+ assert len(ranges) > 1
+ lows = [x for x in ranges[:-1]]
+ highs = [x for x in ranges[1:]]
+ myranges = []
+ for i, pair in enumerate(zip(lows, highs)):
+ if i % 2 == 0:
+ myranges.append(pair)
+ assert len(myranges) == len(ranges) // 2
+ self.ranges = myranges
+
+ def _getRanges(self):
+ return self._ranges
+
+ def _setRanges(self, ranges):
+ # TODO: Check that ranges don't overlap
+ self._ranges = ranges
+ self._totalwidth = sum(r[1] - r[0] for r in ranges)
+
+ runningwidth = 0.0
+ self._divisions = [0.0]
+ for r in ranges:
+ assert(r[1] >= r[0])
+ runningwidth += float(r[1] - r[0])
+ self._divisions.append(runningwidth)
+ self._totalwidth = runningwidth
+ for i in range(len(self._divisions)):
+ self._divisions[i] = float(self._divisions[i]) / float(self._totalwidth)
+
+ ranges = property(_getRanges, _setRanges)
+
+ def _map_unit_to_val(self, x):
+ assert(x >= 0 and x <= 1)
+ idx = None
+ rem = None
+ for i in range(len(self._divisions) - 1):
+ if x >= self._divisions[i] and x < self._divisions[i+1]:
+ idx = i
+ rem = x - self._divisions[i]
+ break
+ if idx is None:
+ raise ValueError("No matching division found in unit interval! How?")
+ val = self.ranges[idx][0] + self._totalwidth * rem
+ return val
+
+ def shoot(self):
+ rand = random.random()
+ val = self._map_unit_to_val(rand)
+ return val
+
+
+class LogSampler(ContinuousSampler):
+ "Randomly sample from an exponential distribution (i.e. uniformly on a log scale)."
+
+ def __init__(self, low, high):
+ self.low = float(low)
+ self.high = float(high)
+
+ def shoot(self):
+ rand = random.random()
+ logval = rand * math.log(self.high) + (1 - rand) * math.log(self.low)
+ val = math.exp(logval)
+ return val
+
+
+class GaussianSampler(ContinuousSampler):
+ "Randomly sample from a 1D Gaussian distribution."
+
+ def __init__(self, mean, sigma):
+ self.mean = float(mean)
+ self.sigma = float(sigma)
+
+ def shoot(self):
+ return random.gauss(self.mean, self.sigma)
+
+
+class InvSampler(ContinuousSampler):
+ "Randomly sample from a 1/x distribution."
+
+ def __init__(self, low, high):
+ self.low = float(low)
+ self.high = float(high)
+
+ def shoot(self):
+ invx = random.uniform(1/self.high, 1/self.low) #< limit inversion not actually necessary
+ return 1./invx
+
+
+########################################
+
+
+class TH1Sampler(ContinuousSampler):
+ "Randomly sample from a 1D ROOT histogram."
+
+ def __init__(self, *args):
+ self.hist = TH1(*args)
+ if self.hist.GetEntries() < 1:
+ raise Exception("Histogram %s is EMPTY! Cannot sample" % self.hist.GetName())
+
+ def shoot(self):
+ return self.hist.GetRandom()
+
+
+########################################
+
+
+## Discrete sequence samplers
+
+class DiscreteSampler(Sampler):
+ "Base class for samplers from lists of discrete values"
+ pass
+
+
+class RandomSeqSampler(DiscreteSampler):
+ "Uniformly random sample from a list of values."
+
+ def __init__(self, *args):
+ if len(args) == 1:
+ self.sequence = args[0]
+ else:
+ self.sequence = args
+
+ def shoot(self):
+ return random.choice(self.sequence)
+# Alias:
+RndmSeq = RandomSeqSampler
+
+
+class CyclicSeqSampler(DiscreteSampler):
+ "Sequentially sample from a list of values, returning to the beginning once exhausted."
+
+ def __init__(self, *args):
+ if len(args) == 1:
+ self.sequence = args[0]
+ else:
+ self.sequence = args
+ self.index = 0
+
+ def shoot(self):
+ self.index = (self.index + 1) % len(self.sequence)
+ return self.sequence[self.index]
+## Alias:
+Sequence = CyclicSeqSampler
+
+
+########################################
+
+
+## Convenience function for sampler-making from Python literals
+
+def mksampler(x):
+ """
+ Automatically cast the provided object to a sampler type. This is used
+ extensively inside the particle and position samplers, so that the user
+ can pass in a primitive type like a number or list and it will be
+ treated as if the more verbose sampler constructors had been called.
+
+ Behaviour:
+ - if x can be called, i.e. x() is valid, we just return x;
+ - a Python list (square brackets) will be converted to a continuous
+ UniformSampler or DisjointUniformSampler;
+ - a Python tuple (round brackets/parentheses) will be treated
+ as a discrete CyclicSeqSampler;
+ - a Python set (curly brackets/braces) will be treated
+ as a discrete RandomSeqSampler;
+ - otherwise a ConstSampler will be created from x, so that x is
+ returned when the sampler is called.
+ """
+ if hasattr(x, "__call__"):
+ return x
+ elif type(x) is list:
+ # NB: disjoint ranges can be given as nested lists, e.g. [(1,2), (4,5)]
+ if len(x) == 2 and type(x[0]) in (int,float) and type(x[1]) in (int,float):
+ #print "MKSAMPLER: Casting %s to UniformSampler" % str(x)
+ return UniformSampler(*x)
+ elif len(x) > 2 or (len(x) > 0 and type(x[0]) not in (int,float)):
+ #print "MKSAMPLER: Casting %s to DisjointUniformSampler" % str(x)
+ return DisjointUniformSampler(x)
+ if len(x) < 2:
+ raise Exception("Supplied list could not be converted to a continuous sampler")
+ elif type(x) is tuple:
+ #print "MKSAMPLER: Casting %s to CyclicSeqSampler" % str(x)
+ return CyclicSeqSampler(*x)
+ elif type(x) is set:
+ #print "MKSAMPLER: Casting %s to RandomSeqSampler" % str(x)
+ return RandomSeqSampler(*x)
+ else:
+ #print "MKSAMPLER: Casting %s to ConstSampler" % str(x)
+ return ConstSampler(x)
+
+
+########################################
+
+
+## Beam-spot (origin vertex) sampling
+
+class PosSampler(Sampler):
+ """
+ Sampler of position 3-vectors, for modelling a beamspot.
+ """
+
+ def __init__(self, x, y, z, t=0):
+ self.x = x
+ self.y = y
+ self.z = z
+ self.t = t
+
+ @property
+ def x(self):
+ "x position sampler"
+ return self._x
+ @x.setter
+ def x(self, x):
+ self._x = mksampler(x)
+
+ @property
+ def y(self):
+ "y position sampler"
+ return self._y
+ @y.setter
+ def y(self, y):
+ self._y = mksampler(y)
+
+ @property
+ def z(self):
+ "z position sampler"
+ return self._z
+ @z.setter
+ def z(self, z):
+ self._z = mksampler(z)
+
+ @property
+ def t(self):
+ "Time sampler"
+ return self._t
+ @t.setter
+ def t(self, t):
+ self._t = mksampler(t)
+
+ def shoot(self):
+ x = self.x()
+ y = self.y()
+ z = self.z()
+ t = self.t()
+ #print "POS =", x, y, z, t
+ return ROOT.TLorentzVector(x, y, z, t)
+
+
+# TODO: Make a 3-Gaussian BeamspotSampler
+
+
+## Momentum sampling
+
+class MomSampler(Sampler):
+ """
+ Base class for four-momentum sampling.
+
+ There are many ways to unambiguously specify four-momenta. Not all are sensible/useful,
+ though. The following are implemented here:
+ * M,px,py,pz
+ * E,M,phi,eta
+ * E,M,phi,y
+ * E,M,phi,theta
+ * pT,M,phi,eta
+ * pT,M,phi,y
+ * pT,M,phi,theta
+
+ Possibly the following (not yet implemented) could be useful: let us know if you
+ need one of these:
+ * E,px,py,pz
+ * pT,E,M,phi
+ """
+ pass
+
+
+class NullMomSampler(MomSampler):
+ "A momentum sampler which just returns a null vector with the given mass."
+
+ def __init__(self, mass=0.0):
+ self.mass = mass
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ def shoot(self):
+ v4 = ROOT.TLorentzVector(0, 0, 0, self.mass)
+ return v4
+
+
+class MXYZSampler(MomSampler):
+ "Create a 4-momentum vector from mass, px, py, pz distributions/samplers."
+
+ def __init__(self, px, py, pz, mass=0.0):
+ self.mass = mass
+ self.px = px
+ self.py = py
+ self.pz = pz
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def px(self):
+ "px sampler"
+ return self._px
+ @px.setter
+ def px(self, x):
+ self._px = mksampler(x)
+
+ @property
+ def py(self):
+ "py sampler"
+ return self._py
+ @py.setter
+ def py(self, x):
+ self._py = mksampler(x)
+
+ @property
+ def pz(self):
+ "pz sampler"
+ return self._pz
+ @pz.setter
+ def pz(self, x):
+ self._pz = mksampler(x)
+
+ def shoot(self):
+ m = self.mass()
+ px = self.px()
+ py = self.py()
+ pz = self.pz()
+ e = math.sqrt(px**2 + py**2 + pz**2 + m**2)
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class EEtaMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from E, eta, m and phi distributions/samplers."
+
+ # TODO: ensure that E >= m!
+
+ def __init__(self, energy, eta, mass=0.0, phi=[0, TWOPI]):
+ self.energy = energy
+ self.eta = eta
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def energy(self):
+ "Energy sampler"
+ return self._e
+ @energy.setter
+ def energy(self, x):
+ self._e = mksampler(x)
+
+ @property
+ def eta(self):
+ "Pseudorapidity sampler"
+ return self._eta
+ @eta.setter
+ def eta(self, x):
+ self._eta = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ eta = - ln(tan(theta/2)) / 2
+ => theta = 2 atan( exp(-eta) )
+ """
+ eta = self.eta()
+ theta = 2 * math.atan(math.exp(-eta))
+ e = self.energy()
+ m = self.mass()
+ p = math.sqrt( e**2 - m**2 )
+ pz = p * math.cos(theta)
+ pt = p * math.sin(theta)
+ phi = self.phi()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class ERapMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from E, y, m and phi distributions."
+
+ # TODO: ensure that E >= m!
+
+ def __init__(self, energy, rap, mass=0.0, phi=[0, TWOPI]):
+ self.energy = energy
+ self.rap = rap
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def energy(self):
+ "Energy sampler"
+ return self._e
+ @energy.setter
+ def energy(self, x):
+ self._e = mksampler(x)
+
+ @property
+ def rap(self):
+ "Rapidity sampler"
+ return self._rap
+ @rap.setter
+ def rap(self, x):
+ self._rap = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ y = 0.5 * ln((E+pz)/(E-pz))
+ -> (E^2 - pz^2) exp(2y) = (E+pz)^2
+ & (E^2 - pz^2) exp(-2y) = (E-pz)^2
+ -> E = sqrt(pt^2 + m^2) cosh(y)
+ -> pz = sqrt(pt^2 + m^2) sinh(y)
+ -> sqrt(pt^2 + m^2) = E / cosh(y)
+ """
+ e = self.energy()
+ y = self.rap()
+ sqrt_pt2_m2 = e / math.cosh(y)
+ pz = sqrt_pt2_m2 * math.sinh(y)
+ m = self.mass()
+ pt = math.sqrt( sqrt_pt2_m2**2 - m**2 )
+ phi = self.phi()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class EThetaMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from E, theta, m and phi distributions/samplers."
+
+ # TODO: ensure that E >= m!
+
+ def __init__(self, energy, theta, mass=0.0, phi=[0, TWOPI]):
+ self.energy = energy
+ self.theta = theta
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def energy(self):
+ "Energy sampler"
+ return self._e
+ @energy.setter
+ def energy(self, x):
+ self._e = mksampler(x)
+
+ @property
+ def theta(self):
+ "Polar angle sampler"
+ return self._theta
+ @theta.setter
+ def theta(self, x):
+ self._theta = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ p = sqrt(e^2 - m^2)
+ pz = p cos(theta)
+ pt = p sin(theta)
+ """
+ e = self.energy()
+ 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()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class PtEtaMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from pt, eta, m and phi distributions/samplers."
+
+ def __init__(self, pt, eta, mass=0.0, phi=[0, TWOPI]):
+ self.pt = pt
+ self.eta = eta
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def pt(self):
+ "Transverse momentum sampler"
+ return self._pt
+ @pt.setter
+ def pt(self, x):
+ self._pt = mksampler(x)
+
+ @property
+ def eta(self):
+ "Pseudorapidity sampler"
+ return self._eta
+ @eta.setter
+ def eta(self, x):
+ self._eta = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ eta = - ln(tan(theta/2)) / 2
+ => theta = 2 atan( exp(-eta) )
+ """
+ eta = self.eta()
+ theta = 2 * math.atan(math.exp(-eta))
+ pt = self.pt()
+ p = pt / math.sin(theta)
+ phi = self.phi()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ pz = p * math.cos(theta)
+ m = self.mass()
+ e = math.sqrt( p**2 + m**2 )
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class PtRapMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from pt, y, m and phi distributions/samplers."
+
+ def __init__(self, pt, rap, mass=0.0, phi=[0, TWOPI]):
+ self.pt = pt
+ self.rap = rap
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def pt(self):
+ "Transverse momentum sampler"
+ return self._pt
+ @pt.setter
+ def pt(self, x):
+ self._pt = mksampler(x)
+
+ @property
+ def rap(self):
+ "Rapidity sampler"
+ return self._rap
+ @rap.setter
+ def rap(self, x):
+ self._rap = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ y = 0.5 * ln((E+pz)/(E-pz))
+ -> (E^2 - pz^2) exp(2y) = (E+pz)^2
+ & (E^2 - pz^2) exp(-2y) = (E-pz)^2
+ -> E = sqrt(pt^2 + m^2) cosh(y)
+ -> pz = sqrt(pt^2 + m^2) sinh(y)
+ -> sqrt(pt^2 + m^2) = E / cosh(y)
+ """
+ pt = self.pt()
+ assert pt >= 0
+ m = self.mass()
+ assert m >= 0
+ sqrt_pt2_m2 = math.sqrt( pt**2 + m**2 )
+ y = self.rap()
+ e = sqrt_pt2_m2 * math.cosh(y)
+ pz = sqrt_pt2_m2 * math.sinh(y)
+ phi = self.phi()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+class PtThetaMPhiSampler(MomSampler):
+ "Create a 4-momentum vector from pt, theta, m and phi distributions/samplers."
+
+ def __init__(self, pt, theta, mass=0.0, phi=[0, TWOPI]):
+ self.pt = pt
+ self.theta = theta
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def pt(self):
+ "Transverse momentum sampler"
+ return self._pt
+ @pt.setter
+ def pt(self, x):
+ self._pt = mksampler(x)
+
+ @property
+ def theta(self):
+ "Polar angle sampler"
+ return self._theta
+ @theta.setter
+ def theta(self, x):
+ self._theta = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ """
+ p = pt / math.sin(theta)
+ pz = p cos(theta)
+ pt = p sin(theta)
+ E = sqrt(p^2 + m^2)
+ """
+ theta = self.theta()
+ pt = self.pt()
+ p = pt / math.sin(theta)
+ phi = self.phi()
+ px = pt * math.cos(phi)
+ py = pt * math.sin(phi)
+ pz = p * math.cos(theta)
+ m = self.mass()
+ e = math.sqrt( p**2 + m**2 )
+ v4 = ROOT.TLorentzVector(px, py, pz, e)
+ return v4
+
+
+# TODO: add the missing ways to specify/sample 4-momenta
+
+
+###########################################################
+
+
+## Combined samplers returning a particle configuration
+
+
+## A default dictionary of particle masses (in MeV)
+MASSES = { 22 : 0.0, # photon
+ 11 : 0.5, # electron
+ 12 : 0.0, # nu_e
+ 13 : 105.7, # muon
+ 14 : 0.0, # nu_mu
+ 15 : 1777.8, # tau
+ 16 : 0.0, # nu_tau
+ 2212 : 938.0, # proton
+ 2112 : 940.0, # neutron
+ 111 : 135.0, # pi0
+ 211 : 140.0, # pi+-
+ 221 : 547.0, # eta
+ 321 : 494.0, # K+-
+ 311 : 598.0 # K0
+ }
+
+
+class SampledParticle(object):
+ """
+ A particle object for use as a return value from the particle samplers.
+ """
+ def __init__(self, pid=None, mom=ROOT.TLorentzVector(0,0,0,0), pos=ROOT.TLorentzVector(0,0,0,0)):
+ """
+ Constructor/initializer: PID is the (int) PDG particle ID code
+ of this particle, mom is its momentum 4-vector, and pos is
+ the vertex 4-position (both as ROOT.TLorentzVector, in MeV).
+ """
+ self.pid = pid
+ self.mom = mom
+ self.pos = pos
+ self.mass = None
+
+
+class ParticleSampler(Sampler):
+ """
+ A simple N-independent-particle sampler.
+ """
+
+ def __init__(self, pid=999,
+ mom=NullMomSampler(),
+ n=1,
+ pos=PosSampler(0, 0, 0)):
+ self.pid = pid
+ self.mom = mom
+ self.n = n
+ self.pos = pos
+ self.massdict = MASSES
+ self.mass_override = True
+
+ @property
+ def pid(self):
+ "Particle ID code sampler"
+ return self._pid
+ @pid.setter
+ def pid(self, x):
+ self._pid = mksampler(x)
+
+ @property
+ def n(self):
+ "Particle number sampler"
+ return self._n
+ @n.setter
+ def n(self, x):
+ self._n = mksampler(x)
+
+ def shoot(self):
+ "Return a vector of sampled particles"
+ numparticles = self.n()
+ rtn = []
+ for i in range(numparticles):
+ ## Sample the particle ID and create a particle
+ pid = self.pid()
+ p = SampledParticle(pid)
+ ## Pass mass info to the v4 sampler and set same generated mass
+ if self.mass_override and abs(pid) in self.massdict:
+ m = self.massdict[abs(pid)]
+ self.mom.mass = m
+ p.mass = m
+ # TODO: Should the particle generated_mass be set from the sampler by default?
+ ## Sample momentum and vertex positions into the particle
+ p.mom = self.mom()
+ p.pos = self.pos()
+ ## Add particle to output list
+ rtn.append(p)
+ return rtn
diff --git a/Generators/ParticleGun/share/common/ParticleGun_Common.py b/Generators/ParticleGun/share/common/ParticleGun_Common.py
new file mode 100644
index 0000000000000000000000000000000000000000..3fab6fb0b7c4093fe92bd0e34e57bc7093b0281a
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_Common.py
@@ -0,0 +1,4 @@
+## Common setup for ParticleGun
+import ParticleGun as PG
+genSeq += PG.ParticleGun()
+evgenConfig.generators += ["ParticleGun"]
diff --git a/Generators/ParticleGun/share/common/ParticleGun_EoverP_Config.py b/Generators/ParticleGun/share/common/ParticleGun_EoverP_Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..8b78a953f31c253a9ead42158bf4d8b0dab77ce0
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_EoverP_Config.py
@@ -0,0 +1,66 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for E/p event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["zach.marshall@cern.ch"]
+
+import ParticleGun as PG
+import ROOT
+from ParticleGun.samplers import *
+class PEtaSampler(PG.MomSampler):
+ "Create a 4-momentum vector from pt, eta, m and phi distributions/samplers."
+
+ def __init__(self, momentum, eta, pid=211, phi=[0, math.pi*2.]):
+ self.momentum = momentum
+ self.eta = eta
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ mass = pdg_table.GetParticle(pid).Mass()
+ self.mass = mass
+ self.phi = phi
+
+ @property
+ def momentum(self):
+ "Momentum sampler"
+ return self._momentum
+ @momentum.setter
+ def momentum(self, x):
+ self._momentum = mksampler(x)
+
+ @property
+ def eta(self):
+ "Pseudorapidity sampler"
+ return self._eta
+ @eta.setter
+ def eta(self, x):
+ self._eta = mksampler(x)
+
+ @property
+ def mass(self):
+ "Mass sampler"
+ return self._m
+ @mass.setter
+ def mass(self, x):
+ self._m = mksampler(x)
+
+ @property
+ def phi(self):
+ "Azimuthal angle sampler"
+ return self._phi
+ @phi.setter
+ def phi(self, x):
+ self._phi = mksampler(x)
+
+ def shoot(self):
+ v4 = ROOT.TLorentzVector()
+ pt = p / math.cosh(self.eta())
+ v4.SetPtEtaPhiM(pt, self.eta(), self.phi(), self.mass())
+ return v4
+
+a_particle = int(jofile.split('_')[-1].split('.py')[0].replace('m','-'))
+
+pg = PG.ParticleGun()
+pg.sampler.pid = int(a_particle) #PID
+pg.sampler.mom = PEtaSampler(momentum=(500,800,1000,1200,1500,2000,3000,4000,5000,6000,7000,8000,9000,10000,11000,13000,15000,17000,20000,\
+ 25000,35000,50000,75000,100000,200000,350000,500000), eta=[-0.3,0.3], pid=int(a_particle))
+genSeq += pg
+
diff --git a/Generators/ParticleGun/share/common/ParticleGun_FastCalo_ChargeFlip_Config.py b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_ChargeFlip_Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..a5399a64019b935aa729244f88c4b5c0ebe5b35f
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_ChargeFlip_Config.py
@@ -0,0 +1,78 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for FastCaloSim event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["david.sosa@cern.ch"]
+
+import ParticleGun as PG
+import ROOT
+
+class MyParticleSampler(PG.ParticleSampler):
+ def __init__(self,energy,eta,pid,shift_z=0):
+ self.pid = pid
+ self.shift_z = shift_z
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ mass = pdg_table.GetParticle(self.pid()).Mass()
+ self.mom1 = PG.EEtaMPhiSampler(energy=energy,eta=eta,mass=mass)
+
+ def shoot(self):
+ pid = self.pid()
+
+ shift_z = self.shift_z
+
+ mom = self.mom1.shoot()
+ pos_temp = mom.Vect().Unit()
+
+ # Would it hit the barrel, or the endcap?
+ if abs(pos_temp.Z())/3550.<pos_temp.Perp()/1148.: # Hit the barrel!
+ pos_temp *= 1148./pos_temp.Perp()
+ else: # Hit the endcap!
+ pos_temp *= 3550./abs(pos_temp.Z())
+
+ # Shift position of vector in the Z direction
+ pos_temp_2 = ROOT.TVector3()
+ pos_temp_2.SetXYZ(pos_temp.X(), pos_temp.Y(), pos_temp.Z()+shift_z)
+ pos_temp_2 *= 1. / pos_temp_2.Mag(); # reduce magnitude of vector
+
+ # recalculate; Would it hit the barrel, or the endcap?
+ if abs(pos_temp_2.Z())/3550.<pos_temp_2.Perp()/1148.:
+ pos_temp_2 *= 1148./pos_temp_2.Perp()
+ else:
+ pos_temp_2 *= 3550./abs(pos_temp_2.Z())
+
+ pos = ROOT.TLorentzVector(pos_temp_2.X(),pos_temp_2.Y(),pos_temp_2.Z(), pos_temp_2.Mag())
+
+ #print "pid ",pid
+
+ return [ PG.SampledParticle( pid , mom , pos ) ]
+
+myE = float(jofile.split('_E')[1].split('_')[0])
+myZV = float(jofile.split('_')[-1].split('.py')[0].replace("m","-"))
+
+myPDGID = jofile.split('_pid')[1].split('_')[0].replace('n','-')
+myPDGID = int(float(myPDGID.replace('p','')))
+
+eta_li = []
+
+if "disj" in jofile:
+ myLowEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[0].replace('m','-'))
+ myLowEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[1].replace('m','-'))
+ myHighEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[2].replace('m','-'))
+ myHighEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[3].replace('m','-'))
+ eta_li.extend([myLowEta1,myLowEta2,myHighEta1,myHighEta2])
+
+else:
+ myLowEta = 0.01*float(jofile.split('eta')[1].split('_')[0].replace('m','-'))
+ myHighEta = 0.01*float(jofile.split('eta')[1].split('_')[1].replace('m','-'))
+ eta_li.extend([myLowEta,myHighEta])
+
+
+print "================ SETTTINGS ================="
+print ("energy = ", myE)
+print ("eta = ", eta_li)
+print ("pid = ", myPDGID)
+print ("shift_z = ", myZV)
+print "============================================"
+
+genSeq += PG.ParticleGun()
+genSeq.ParticleGun.sampler = MyParticleSampler(energy=myE,eta=eta_li,pid=(myPDGID,myPDGID),shift_z=myZV) #unmixed
diff --git a/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config.py b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..1b2e9a68bc5e1c1612bf2e294c58dcc472f700fb
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config.py
@@ -0,0 +1,101 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for FastCaloSim event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["david.sosa@cern.ch"]
+
+import ParticleGun as PG
+import ROOT
+
+class MyParticleSampler(PG.ParticleSampler):
+ def __init__(self,energy,eta,pid,shift_z=0):
+ self.pid = pid
+ self.shift_z = shift_z
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ mass = pdg_table.GetParticle(self.pid()).Mass()
+ self.mom1 = PG.EEtaMPhiSampler(energy=energy,eta=eta,mass=mass)
+
+ def shoot(self):
+ pid = self.pid()
+
+ shift_z = self.shift_z
+
+ mom = self.mom1.shoot()
+ pos_temp = mom.Vect().Unit()
+
+ # Define geometry
+ barrelR1 = 1148.0
+ barrelR2 = 120.0
+ barrelR3 = 41.0
+ endcapZ1 = 3550.0
+ endcapZ2 = 4587.0
+
+ # Would it hit the barrel, or the endcap?
+ tanTheta = pos_temp.Perp() / abs( pos_temp.Z() );
+ if tanTheta > barrelR1 / endcapZ1:
+ pos_temp *= barrelR1 / pos_temp.Perp()
+ elif tanTheta > barrelR2 / endcapZ1:
+ pos_temp *= endcapZ1 / abs( pos_temp.Z() )
+ elif tanTheta > barrelR2 / endcapZ2:
+ pos_temp *= barrelR2 / pos_temp.Perp()
+ elif tanTheta > barrelR3 / endcapZ2:
+ pos_temp *= endcapZ2 / abs( pos_temp.Z() )
+ else:
+ pos_temp *= barrelR3 / pos_temp.Perp()
+
+ # Shift position of vector in the Z direction
+ pos_temp_2 = ROOT.TVector3()
+ pos_temp_2.SetXYZ(pos_temp.X(), pos_temp.Y(), pos_temp.Z()+shift_z)
+ pos_temp_2 *= 1. / pos_temp_2.Mag(); # reduce magnitude of vector
+
+ # recalculate; Would it hit the barrel, or the endcap?
+ tanTheta_2 = pos_temp_2.Perp() / abs( pos_temp_2.Z() );
+ if tanTheta_2 > barrelR1 / endcapZ1:
+ pos_temp_2 *= barrelR1 / pos_temp_2.Perp()
+ elif tanTheta_2 > barrelR2 / endcapZ1:
+ pos_temp_2 *= endcapZ1 / abs( pos_temp_2.Z() )
+ elif tanTheta_2 > barrelR2 / endcapZ2:
+ pos_temp_2 *= barrelR2 / pos_temp_2.Perp()
+ elif tanTheta_2 > barrelR3 / endcapZ2:
+ pos_temp_2 *= endcapZ2 / abs( pos_temp_2.Z() )
+ else:
+ pos_temp_2 *= barrelR3 / pos_temp_2.Perp()
+
+ pos = ROOT.TLorentzVector(pos_temp_2.X(),pos_temp_2.Y(),pos_temp_2.Z(), pos_temp_2.Mag())
+
+ #print "pid ",pid
+
+ return [ PG.SampledParticle( pid , mom , pos ) ]
+
+myE = float(jofile.split('_E')[1].split('_')[0])
+myZV = float(jofile.split('_')[-1].split('.py')[0].replace("m","-"))
+myPDGID = int(float(jofile.split('_pid')[1].split('_')[0].replace('m','-')))
+
+eta_li = []
+
+if "disj" in jofile:
+ myLowEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[0].replace('m','-'))
+ myLowEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[1].replace('m','-'))
+ myHighEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[2].replace('m','-'))
+ myHighEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[3].replace('m','-'))
+ eta_li.extend([myLowEta1,myLowEta2,myHighEta1,myHighEta2])
+
+else:
+ myLowEta = 0.01*float(jofile.split('eta')[1].split('_')[0].replace('m','-'))
+ myHighEta = 0.01*float(jofile.split('eta')[1].split('_')[1].replace('m','-'))
+ eta_li.extend([myLowEta,myHighEta])
+
+
+print "================ SETTTINGS ================="
+print ("energy = ", myE)
+print ("eta = ", eta_li)
+print ("pid = ", myPDGID)
+print ("shift_z = ", myZV)
+print "============================================"
+
+genSeq += PG.ParticleGun()
+if myPDGID != 22:
+ genSeq.ParticleGun.sampler = MyParticleSampler(energy=myE,eta=eta_li,pid=(-myPDGID,myPDGID),shift_z=myZV)
+else:
+ genSeq.ParticleGun.sampler = MyParticleSampler(energy=myE,eta=eta_li,pid=myPDGID,shift_z=myZV)
+
diff --git a/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config_Erange.py b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config_Erange.py
new file mode 100644
index 0000000000000000000000000000000000000000..75ebc0621e7ba4e6803ae25c20a6d7438fc45466
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_Config_Erange.py
@@ -0,0 +1,103 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for FastCaloSim event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["david.sosa@cern.ch"]
+
+import ParticleGun as PG
+import ROOT
+
+class MyParticleSampler(PG.ParticleSampler):
+ def __init__(self,energy,eta,pid,shift_z=0):
+ self.pid = pid
+ self.shift_z = shift_z
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ mass = pdg_table.GetParticle(self.pid()).Mass()
+ self.mom1 = PG.EEtaMPhiSampler(energy=energy,eta=eta,mass=mass)
+
+ def shoot(self):
+ pid = self.pid()
+
+ shift_z = self.shift_z
+
+ mom = self.mom1.shoot()
+ pos_temp = mom.Vect().Unit()
+
+ # Define geometry
+ barrelR1 = 1148.0
+ barrelR2 = 120.0
+ barrelR3 = 41.0
+ endcapZ1 = 3550.0
+ endcapZ2 = 4587.0
+
+ # Would it hit the barrel, or the endcap?
+ tanTheta = pos_temp.Perp() / abs( pos_temp.Z() );
+ if tanTheta > barrelR1 / endcapZ1:
+ pos_temp *= barrelR1 / pos_temp.Perp()
+ elif tanTheta > barrelR2 / endcapZ1:
+ pos_temp *= endcapZ1 / abs( pos_temp.Z() )
+ elif tanTheta > barrelR2 / endcapZ2:
+ pos_temp *= barrelR2 / pos_temp.Perp()
+ elif tanTheta > barrelR3 / endcapZ2:
+ pos_temp *= endcapZ2 / abs( pos_temp.Z() )
+ else:
+ pos_temp *= barrelR3 / pos_temp.Perp()
+
+ # Shift position of vector in the Z direction
+ pos_temp_2 = ROOT.TVector3()
+ pos_temp_2.SetXYZ(pos_temp.X(), pos_temp.Y(), pos_temp.Z()+shift_z)
+ pos_temp_2 *= 1. / pos_temp_2.Mag(); # reduce magnitude of vector
+
+ # recalculate; Would it hit the barrel, or the endcap?
+ tanTheta_2 = pos_temp_2.Perp() / abs( pos_temp_2.Z() );
+ if tanTheta_2 > barrelR1 / endcapZ1:
+ pos_temp_2 *= barrelR1 / pos_temp_2.Perp()
+ elif tanTheta_2 > barrelR2 / endcapZ1:
+ pos_temp_2 *= endcapZ1 / abs( pos_temp_2.Z() )
+ elif tanTheta_2 > barrelR2 / endcapZ2:
+ pos_temp_2 *= barrelR2 / pos_temp_2.Perp()
+ elif tanTheta_2 > barrelR3 / endcapZ2:
+ pos_temp_2 *= endcapZ2 / abs( pos_temp_2.Z() )
+ else:
+ pos_temp_2 *= barrelR3 / pos_temp_2.Perp()
+
+ pos = ROOT.TLorentzVector(pos_temp_2.X(),pos_temp_2.Y(),pos_temp_2.Z(), pos_temp_2.Mag())
+
+ #print "pid ",pid
+
+ return [ PG.SampledParticle( pid , mom , pos ) ]
+
+E_li = []
+myLowE = float(jofile.split('_E')[1].split('_')[0])
+myHighE = float(jofile.split('_E')[1].split('_')[1])
+E_li.extend([myLowE,myHighE])
+
+myZV = float(jofile.split('_')[-1].split('.py')[0].replace("m","-"))
+myPDGID = int(float(jofile.split('_pid')[1].split('_')[0].replace('m','-')))
+
+eta_li = []
+
+if "disj" in jofile:
+ myLowEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[0].replace('m','-'))
+ myLowEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[1].replace('m','-'))
+ myHighEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[2].replace('m','-'))
+ myHighEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[3].replace('m','-'))
+ eta_li.extend([myLowEta1,myLowEta2,myHighEta1,myHighEta2])
+
+else:
+ myLowEta = 0.01*float(jofile.split('eta')[1].split('_')[0].replace('m','-'))
+ myHighEta = 0.01*float(jofile.split('eta')[1].split('_')[1].replace('m','-'))
+ eta_li.extend([myLowEta,myHighEta])
+
+
+print "================ SETTTINGS ================="
+print ("energy = ", E_li)
+print ("eta = ", eta_li)
+print ("pid = ", myPDGID)
+print ("shift_z = ", myZV)
+print "============================================"
+
+genSeq += PG.ParticleGun()
+print "E_li = ", E_li, ", eta_li = ", eta_li, ", pid = ", myPDGID, ", myZV = ", myZV
+genSeq.ParticleGun.sampler = MyParticleSampler(energy=E_li,eta=eta_li,pid=myPDGID,shift_z=myZV)
+
diff --git a/Generators/ParticleGun/share/common/ParticleGun_FastCalo_NoChargeFlip_Config.py b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_NoChargeFlip_Config.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ba60ef2bc9d3a9195ada72e9a2232864f173567
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_FastCalo_NoChargeFlip_Config.py
@@ -0,0 +1,78 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for FastCaloSim event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["david.sosa@cern.ch"]
+
+import ParticleGun as PG
+import ROOT
+
+class MyParticleSampler(PG.ParticleSampler):
+ def __init__(self,energy,eta,pid,shift_z=0):
+ self.pid = pid
+ self.shift_z = shift_z
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ mass = pdg_table.GetParticle(self.pid()).Mass()
+ self.mom1 = PG.EEtaMPhiSampler(energy=energy,eta=eta,mass=mass)
+
+ def shoot(self):
+ pid = self.pid()
+
+ shift_z = self.shift_z
+
+ mom = self.mom1.shoot()
+ pos_temp = mom.Vect().Unit()
+
+ # Would it hit the barrel, or the endcap?
+ if abs(pos_temp.Z())/3550.<pos_temp.Perp()/1148.: # Hit the barrel!
+ pos_temp *= 1148./pos_temp.Perp()
+ else: # Hit the endcap!
+ pos_temp *= 3550./abs(pos_temp.Z())
+
+ # Shift position of vector in the Z direction
+ pos_temp_2 = ROOT.TVector3()
+ pos_temp_2.SetXYZ(pos_temp.X(), pos_temp.Y(), pos_temp.Z()+shift_z)
+ pos_temp_2 *= 1. / pos_temp_2.Mag(); # reduce magnitude of vector
+
+ # recalculate; Would it hit the barrel, or the endcap?
+ if abs(pos_temp_2.Z())/3550.<pos_temp_2.Perp()/1148.:
+ pos_temp_2 *= 1148./pos_temp_2.Perp()
+ else:
+ pos_temp_2 *= 3550./abs(pos_temp_2.Z())
+
+ pos = ROOT.TLorentzVector(pos_temp_2.X(),pos_temp_2.Y(),pos_temp_2.Z(), pos_temp_2.Mag())
+
+ #print "pid ",pid
+
+ return [ PG.SampledParticle( pid , mom , pos ) ]
+
+myE = float(jofile.split('_E')[1].split('_')[0])
+myZV = float(jofile.split('_')[-1].split('.py')[0].replace("m","-"))
+
+myPDGID = jofile.split('_pid')[1].split('_')[0].replace('n','-')
+myPDGID = int(float(myPDGID.split('_pid')[1].split('_')[0].replace('p','')))
+
+eta_li = []
+
+if "disj" in jofile:
+ myLowEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[0].replace('m','-'))
+ myLowEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[1].replace('m','-'))
+ myHighEta1 = 0.01*float(jofile.split('eta_')[1].split('_')[2].replace('m','-'))
+ myHighEta2 = 0.01*float(jofile.split('eta_')[1].split('_')[3].replace('m','-'))
+ eta_li.extend([myLowEta1,myLowEta2,myHighEta1,myHighEta2])
+
+else:
+ myLowEta = 0.01*float(jofile.split('eta')[1].split('_')[0].replace('m','-'))
+ myHighEta = 0.01*float(jofile.split('eta')[1].split('_')[1].replace('m','-'))
+ eta_li.extend([myLowEta,myHighEta])
+
+
+print "================ SETTTINGS ================="
+print ("energy = ", myE)
+print ("eta = ", eta_li)
+print ("pid = ", myPDGID)
+print ("shift_z = ", myZV)
+print "============================================"
+
+genSeq += PG.ParticleGun()
+genSeq.ParticleGun.sampler = MyParticleSampler(energy=myE,eta=eta_li,pid=(myPDGID,myPDGID),shift_z=myZV) #unmixed
diff --git a/Generators/ParticleGun/share/common/ParticleGun_SamplingFraction.py b/Generators/ParticleGun/share/common/ParticleGun_SamplingFraction.py
new file mode 100644
index 0000000000000000000000000000000000000000..54557f0d5d2afacedb09999acac22e02ad8576ac
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_SamplingFraction.py
@@ -0,0 +1,97 @@
+#! -*- python -*-
+evgenConfig.description = "Single particle gun for Sampling Fraction event generation"
+evgenConfig.keywords = ["singleParticle",]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["michael.duehrssen@cern.ch"]
+
+import ParticleGun as PG
+import ROOT, math, random
+
+class MyParticleSampler(PG.ParticleSampler):
+ """
+ Projective showers starting at entrance of calorimeter, flat in eta, constant energy
+ """
+
+ def __init__(self,pid=11,momentum=50000.,eta1=0.,eta2=1.4,bec=0,radius=1500.,z=3740.5):
+ self.pid = pid
+ self.momentum = momentum
+ self.eta1 = eta1
+ self.eta2 = eta2
+ pdg_table = ROOT.TDatabasePDG.Instance()
+ self.mass = pdg_table.GetParticle(self.pid()).Mass()
+ self.bec=bec
+ self.radius=radius
+ self.z=z
+
+ def shoot(self):
+ rtn=[]
+ eta = random.uniform(self.eta1, self.eta2)
+ phi = random.uniform(0, math.tau) # tau = 2 * pi
+ v4 = ROOT.TLorentzVector()
+ pt = self.momentum / math.cosh(eta)
+ v4.SetPtEtaPhiM(pt, eta, phi, self.mass)
+ if self.bec==0:
+ radius= self.radius
+ x=radius*math.cos(phi)
+ y=radius*math.sin(phi)
+ z=radius*math.sinh(eta)
+ else:
+ z=self.z
+ radius=z/math.sinh(eta)
+ x=radius*math.cos(phi)
+ y=radius*math.sin(phi)
+ t=math.sqrt(x*x+y*y+z*z)
+ vp = ROOT.TLorentzVector(x,y,z,t)
+ p = PG.SampledParticle(pid=self.pid(),mom=v4,pos=vp)
+ #print "E,eta,phi,mass ",e,eta,phi,self.mass," position ",x,y,z," pid=",p.pid
+ rtn.append(p)
+ return rtn
+
+##MC15 style with Generate_tf.py
+#args=jofile.split('.py')[0]
+
+##MC16 style with Gen_tf.py
+FIRST_DIR = (os.environ['JOBOPTSEARCHPATH']).split(":")[0]
+jofiles = [f for f in os.listdir(FIRST_DIR) if (f.startswith('mc') and f.endswith('.py'))]
+
+print "================ SETTTINGS ================="
+print ("jofiles = ", jofiles)
+
+### parse options from MC job-options filename
+args = jofiles[0].split('.py')[0]
+print ("args = ", args)
+
+myMomentum = float(args.split('_Mom')[1].split('_')[0])
+print ("Momentum = ", myMomentum,"MeV")
+
+myPDGID = int(float(args.split('_pid')[1].split('_')[0].replace('m','-')))
+print ("pid = ", myPDGID)
+
+myLowEta = 0.01*float(args.split('eta_')[1].split('_')[0].replace('m','-'))
+print ("etalow = ", myLowEta)
+
+myHighEta = 0.01*float(args.split('eta_')[1].split('_')[1].replace('m','-'))
+print ("etahigh = ", myHighEta)
+
+if "_Radius" in args:
+ myRadius = 0.001*float(args.split('_Radius')[1].split('_')[0]) #Argument needs to by in mum, since a "." in the filename is not allowed
+else:
+ myRadius = 1500.
+print ("radius = ", myRadius,"mm")
+
+if "_Z" in args:
+ myZ = 0.001*float(args.split('_Z')[1].split('_')[0]) #Argument needs to by in mum, since a "." in the filename is not allowed
+else:
+ myZ = 3740.5
+print ("Z = ", myZ,"mm")
+
+if "bec" in args:
+ bec=1
+else:
+ bec=0
+print ("bec = ", bec)
+print "============================================"
+
+genSeq += PG.ParticleGun()
+genSeq.ParticleGun.sampler = MyParticleSampler(momentum=myMomentum,eta1=myLowEta,eta2=myHighEta,pid=myPDGID,bec=bec,radius=myRadius,z=myZ)
+
diff --git a/Generators/ParticleGun/share/common/ParticleGun_SingleHECO.py b/Generators/ParticleGun/share/common/ParticleGun_SingleHECO.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ffb36dd28235948f0ec1d298c121902267f27d1
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_SingleHECO.py
@@ -0,0 +1,74 @@
+
+PDG = 10000000 + int(float(charge)*100.0)
+loE = (float(mass) + 10.)*1000.
+hiE = (float(mass) + 6000.)*1000.
+MeVmass=float(mass)*1000.
+#--------------------------------------------------------------
+# Configuration for EvgenJobTransforms
+#--------------------------------------------------------------
+evgenConfig.description = "Single HECO generation for Mass=%s, Charge=%s in MC15" % (mass,charge)
+evgenConfig.keywords = ["exotic", "singleParticle","highElectricChargeObject"]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["anlionti@cern.ch"]
+
+evgenConfig.specialConfig = 'MASS=%s;CHARGE=%s;preInclude=SimulationJobOptions/preInclude.Qball.py' % (mass,charge)
+
+
+
+#--------------------------------------------------------------
+# Configuration for ParticleGun
+#--------------------------------------------------------------
+include("ParticleGun/ParticleGun_Common.py")
+
+import ParticleGun as PG
+PG.MASSES[PDG] = float(MeVmass)
+genSeq.ParticleGun.sampler.pid = (-PDG, PDG)
+genSeq.ParticleGun.sampler.mom = PG.EEtaMPhiSampler(energy=[loE,hiE], eta=[-2,2])
+
+
+#--------------------------------------------------------------
+# Edit PDGTABLE.MeV with monopole mass
+#--------------------------------------------------------------
+ALINE1="M %s %s.E+03 +0.0E+00 -0.0E+00 Monopole 0" % (PDG,mass)
+ALINE2="W %s 0.E+00 +0.0E+00 -0.0E+00 Monopole 0" % (PDG)
+
+import os
+import sys
+
+pdgmod = os.path.isfile('PDGTABLE.MeV')
+if pdgmod is True:
+ os.remove('PDGTABLE.MeV')
+os.system('get_files -data PDGTABLE.MeV')
+f=open('PDGTABLE.MeV','a')
+f.writelines(str(ALINE1))
+f.writelines('\n')
+f.writelines(str(ALINE2))
+f.writelines('\n')
+f.close()
+
+del ALINE1
+del ALINE2
+
+#--------------------------------------------------------------
+# Edit G4particle_whitelist.txt with monopole
+#--------------------------------------------------------------
+
+ALINE1="%s qb %s.E+03 (Mev/c) lepton %s" % (PDG,mass,charge)
+ALINE2="-%s qbbar %s.E+03 (Mev/c) lepton -%s" % (PDG,mass,charge)
+
+import os
+import sys
+
+pdgmod = os.path.isfile('G4particle_whitelist.txt')
+if pdgmod is True:
+ os.remove('G4particle_whitelist.txt')
+os.system('get_files -data G4particle_whitelist.txt')
+f=open('G4particle_whitelist.txt','a')
+f.writelines(str(ALINE1))
+f.writelines('\n')
+f.writelines(str(ALINE2))
+f.writelines('\n')
+f.close()
+
+del ALINE1
+del ALINE2
diff --git a/Generators/ParticleGun/share/common/ParticleGun_SingleMonopole.py b/Generators/ParticleGun/share/common/ParticleGun_SingleMonopole.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a6bf0574561af0f4cab4a20e68248dcc3111bae
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_SingleMonopole.py
@@ -0,0 +1,74 @@
+
+PDG = 4110000
+loE = (float(monmass) + 10.)*1000.
+hiE = (float(monmass) + 6000.)*1000.
+MeVmass=float(monmass)*1000.
+#--------------------------------------------------------------
+# Configuration for EvgenJobTransforms
+#--------------------------------------------------------------
+evgenConfig.description = "Single magnetic monopole generation for Mass=%s, Gcharge=%s in MC15" % (monmass,gcharge)
+evgenConfig.keywords = ["exotic", "magneticMonopole", "singleParticle"]
+evgenConfig.generators = ["ParticleGun"]
+evgenConfig.contact = ["anlionti@cern.ch"]
+
+evgenConfig.specialConfig = 'MASS=%s;GCHARGE=%s;preInclude=SimulationJobOptions/preInclude.Monopole.py' % (monmass,gcharge)
+
+
+
+#--------------------------------------------------------------
+# Configuration for ParticleGun
+#--------------------------------------------------------------
+include("ParticleGun/ParticleGun_Common.py")
+
+import ParticleGun as PG
+PG.MASSES[4110000] = float(MeVmass)
+genSeq.ParticleGun.sampler.pid = (-PDG, PDG)
+genSeq.ParticleGun.sampler.mom = PG.EEtaMPhiSampler(energy=[loE,hiE], eta=[-2,2])
+
+
+#--------------------------------------------------------------
+# Edit PDGTABLE.MeV with monopole mass
+#--------------------------------------------------------------
+ALINE1="M 4110000 %s.E+03 +0.0E+00 -0.0E+00 Monopole 0" % (monmass)
+ALINE2="W 4110000 0.E+00 +0.0E+00 -0.0E+00 Monopole 0"
+
+import os
+import sys
+
+pdgmod = os.path.isfile('PDGTABLE.MeV')
+if pdgmod is True:
+ os.remove('PDGTABLE.MeV')
+os.system('get_files -data PDGTABLE.MeV')
+f=open('PDGTABLE.MeV','a')
+f.writelines(str(ALINE1))
+f.writelines('\n')
+f.writelines(str(ALINE2))
+f.writelines('\n')
+f.close()
+
+del ALINE1
+del ALINE2
+
+#--------------------------------------------------------------
+# Edit G4particle_whitelist.txt with monopole
+#--------------------------------------------------------------
+
+ALINE1="4110000 mm %s.E+03 (Mev/c) lepton %s" % (monmass,gcharge)
+ALINE2="-4110000 mmbar %s.E+03 (Mev/c) lepton -%s" % (monmass,gcharge)
+
+import os
+import sys
+
+pdgmod = os.path.isfile('G4particle_whitelist.txt')
+if pdgmod is True:
+ os.remove('G4particle_whitelist.txt')
+os.system('get_files -data G4particle_whitelist.txt')
+f=open('G4particle_whitelist.txt','a')
+f.writelines(str(ALINE1))
+f.writelines('\n')
+f.writelines(str(ALINE2))
+f.writelines('\n')
+f.close()
+
+del ALINE1
+del ALINE2
diff --git a/Generators/ParticleGun/share/common/ParticleGun_egammaET.py b/Generators/ParticleGun/share/common/ParticleGun_egammaET.py
new file mode 100644
index 0000000000000000000000000000000000000000..04ea9b92a3e209afc279d743307dff393f1d8839
--- /dev/null
+++ b/Generators/ParticleGun/share/common/ParticleGun_egammaET.py
@@ -0,0 +1,51 @@
+__doc__ = "Holds a 4-momentum sampler according to the egamma Et spectrum"
+
+import ParticleGun as PG
+from GaudiKernel.SystemOfUnits import GeV
+
+def dbnFermiDirac(x,mu,kT):
+ import math
+ arg = (x-mu)/kT
+ if arg < -20 : # avoid numerical underflows
+ result = 1
+ elif arg > 20 : # avoid numerical overflows
+ result = 0
+ else :
+ div = math.exp(arg)+1
+ result = 1/div
+ return result
+
+class egammaETSampler(PG.PtEtaMPhiSampler):
+ "4-momentum sampler according to the egamma Et spectrum."
+ def __init__(self, pid, eta=[-2.5, 2.5], phi=[0, PG.TWOPI],
+ mu1 = 0.5, kT1 = 0.1, mu2 = 200, kT2 = 20, y0 = 0.005, PtMin = 0 , PtMax = 3e3, nBins=None):
+ """
+ Parameters for the MVA-shaped spectrum : higher density in the < 100 GeV range
+ PtMin = 0 # minimum Pt
+ PtMax = 3000 # maximum Pt (3 TeV)
+ nBins # number of bins (one every 100 MeV by default)
+ mu1 = 0.5 # mu1,kT1 : smooth but steep ramp-up from 0 to 1 GeV (requested by TauCP)
+ kT1 = 0.1
+ mu2 = 200 # mu2,kT2 : smooth, slow ramp-down in the 100-300 GeV range
+ kT2 = 20
+ y0 = 0.005 # y0 : baseline for low-density at high ET up to PtMax
+ """
+ self.m = PG.MASSES[abs(pid)]
+ self.eta = eta
+ self.phi = phi
+
+ # Create and fill a very fine-grained histogram
+ from ROOT import TH1D
+ etSpectrumFullRange = TH1D("ETSpectrumFullRange",
+ "Reference ET spectrum for egamma MVA calib",
+ int(nBins or (PtMax - PtMin)*10), PtMin , PtMax)
+ for i in xrange(etSpectrumFullRange.GetNbinsX()):
+ x = etSpectrumFullRange.GetBinCenter(i+1)
+ y1 = dbnFermiDirac(x,mu1,kT1)
+ y2 = dbnFermiDirac(x,mu2,kT2)
+ y = y0 - y1 + y2
+ etSpectrumFullRange.SetBinContent(i+1,y)
+ self.hist = PG.TH1(etSpectrumFullRange) #< wrap *after* populating
+
+ def pt(self):
+ return self.hist.GetRandom() * GeV
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_constenergy_flateta.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_constenergy_flateta.py
new file mode 100644
index 0000000000000000000000000000000000000000..3c3cb0321b37e7f3d5ec3fcfc2602b2f1cba2cb0
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_constenergy_flateta.py
@@ -0,0 +1,17 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+pg = PG.ParticleGun()
+pg.randomSeed = 123456
+pg.sampler.pid = {11,-11,211,111}
+pg.sampler.mom = PG.EEtaMPhiSampler(energy=10000, eta=[-2,2])
+topSeq += pg
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_correlated.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_correlated.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb0a9437388d3b21fffd744c67cfe590dbddaf4e
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_correlated.py
@@ -0,0 +1,34 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+
+class MyParticleSampler(PG.ParticleSampler):
+ "A special sampler with two _correlated_ particles."
+
+ def __init__(self):
+ self.mom1 = PG.PtEtaMPhiSampler(pt=25000, eta=[-2,2])
+
+ def shoot(self):
+ "Return a vector of sampled particles"
+ p1 = PG.SampledParticle(11, self.mom1.shoot())
+ eta1 = p1.mom.Eta()
+ phi1 = p1.mom.Phi()
+ # TODO: will phi be properly wrapped into range?
+ mom2 = PG.PtEtaMPhiSampler(pt=25000,
+ eta=[eta1-0.5, eta1+0.5],
+ phi=[phi1-0.5, phi1+0.5])
+ p2 = PG.SampledParticle(11, mom2.shoot())
+ return [p1, p2]
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.randomSeed = 123456
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_corrhist.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_corrhist.py
new file mode 100644
index 0000000000000000000000000000000000000000..7c0cd41b88e35f0290b137e5eebe80296a90b3ca
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_corrhist.py
@@ -0,0 +1,41 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+## ROOT 2D histogram sampling alg (in ParticleGun.histsampling) by Andy Buckley
+## Thanks to Alejandro Alonso for the initial Athena example on which this is based.
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+
+class PtEtaHistParticleSampler(PG.ParticleSampler):
+ "Particle sampler with correlated pT and eta from a 2D histogram."
+
+ def __init__(self, pid, histfile, num=100):
+ self.pid = PG.mksampler(pid)
+ self.hist = PG.TH2(histfile, "h_pt_eta")
+ self.numparticles = num
+
+ def shoot(self):
+ "Return a vector of sampled particles from the provided pT--eta histogram"
+ particles = []
+ for i in xrange(self.numparticles):
+ ptrand, etarand = self.hist.GetRandom()
+ ptrand *= 1000 # NB. This _particular_ histogram is in GeV, but Athena needs MeV!
+ # TODO: Provide 4-mom construction functions to avoid building this one-time sampler
+ pid = self.pid()
+ mom = PG.PtEtaMPhiSampler(pt=ptrand, eta=etarand, mass=PG.MASSES[abs(pid)])
+ p = PG.SampledParticle(pid, mom())
+ #print p.mom.Pt(), "\t", p.mom.Eta(), "\t", p.mom.Phi(), "\t", p.mom.M()
+ particles.append(p)
+ return particles
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.randomSeed = 123456
+topSeq.ParticleGun.sampler = PtEtaHistParticleSampler(11, "data_histos_el_1470pt.root")
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatcurvature_flatip.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatcurvature_flatip.py
new file mode 100644
index 0000000000000000000000000000000000000000..38233563de9673b1d9ebd1ae434b57f6051ae70e
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatcurvature_flatip.py
@@ -0,0 +1,41 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+
+class MyParticleSampler(PG.ParticleSampler):
+ """
+ A special sampler to generate single particles flat in 1/pT and in
+ impact parameter to the beam, with flat z0.
+ """
+
+ def __init__(self):
+ psamp = PG.PtEtaMPhiSampler(pt=PG.InvSampler(4000, 400000), eta=[0.1,0.3], phi=[0.3, 0.5])
+ xsamp = PG.PosSampler(0, 0, [-150,150], 0)
+ PG.ParticleSampler.__init__(self, pid={13,-13}, mom=psamp, pos=xsamp)
+ self.ip = PG.mksampler([-2,2])
+
+ def shoot(self):
+ "Return a vector of sampled particles"
+ ps = PG.ParticleSampler.shoot(self)
+ assert len(ps) == 1
+ p = ps[0]
+ from math import sqrt
+ m = -p.mom.X() / p.mom.Y() #< gradient of azimuthal IP sampling line, perp to mom
+ x = self.ip() / sqrt(1 + m**2) #< just decomposing sampled IP into x component...
+ y = m*x #< ... and y-component
+ p.pos.SetX(x)
+ p.pos.SetY(m*x)
+ return [p]
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.randomSeed = 123456
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatpt_2particle.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatpt_2particle.py
new file mode 100644
index 0000000000000000000000000000000000000000..97ed64f8857e82dbcaeb06085fc84d84eb246c0b
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_flatpt_2particle.py
@@ -0,0 +1,20 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+pg = PG.ParticleGun()
+pg.randomSeed = 123456
+pg.samplers.append(PG.ParticleSampler()) # add a second sampler
+pg.samplers[0].pid = (-13, 13) # cycle mu+-
+pg.samplers[0].mom = PG.PtEtaMPhiSampler(pt=[4000, 100000], eta=[1.0, 3.2]) # flat in pt and +ve eta
+pg.samplers[1].pid = (13, -13) # cycle mu-+
+pg.samplers[1].mom = PG.PtEtaMPhiSampler(pt=[4000, 100000], eta=[-3.2, -1.0]) # flat in pt and -ve eta
+topSeq += pg
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_fwd_sequence.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_fwd_sequence.py
new file mode 100644
index 0000000000000000000000000000000000000000..d1d4746fb7b7002fb8e096058c1b2b5d1c4eb348
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_fwd_sequence.py
@@ -0,0 +1,19 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+pg = PG.ParticleGun()
+pg.randomSeed = 123456
+pg.sampler.pid = (2112, 22, 2112, 22)
+pg.sampler.mom = PG.EThetaMPhiSampler(energy=(1360000, 500000, 1360000, 500000),
+ theta=(0, 0, PG.PI, PG.PI))
+pg.sampler.pos = PG.PosSampler(x=[-120,-100], y=[-10,10], z=203950)
+topSeq += pg
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/ParticleGun/share/examples/jobOption.ParticleGun_vtx.py b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_vtx.py
new file mode 100644
index 0000000000000000000000000000000000000000..cdacb7ff5a678cb95894a39a14f0ead507c41a4b
--- /dev/null
+++ b/Generators/ParticleGun/share/examples/jobOption.ParticleGun_vtx.py
@@ -0,0 +1,18 @@
+#! -*- python -*-
+
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 100
+
+import ParticleGun as PG
+pg = PG.ParticleGun()
+pg.randomSeed = 123456
+pg.sampler.pid = 13
+pg.sampler.pos = PG.PosSampler(x=3140.0, y=[-154.134,154.134], z=[4938.76,5121.29], t=5929.7)
+pg.sampler.mom = PG.EEtaMPhiSampler(energy=100000, eta=1.25, phi=0.0)
+topSeq += pg
+
+include("GeneratorUtils/postJO.CopyWeights.py")
+include("GeneratorUtils/postJO.PoolOutput.py")
+include("GeneratorUtils/postJO.DumpMC.py")