diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 8170bbd33b732f1972d3fbdedea626c61948d95d..58f34f08b6cd3d23322e55f7737aa43559a0047f 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -101,9 +101,9 @@ class EvgenValidation(EvgenAnalysisAlg):
self.hists.add("P_M", 1000, 0, 10000)
self.hists.add("Pz_M", 1000, 0, 10000)
self.hists.add("Pt_M", 1000, 0, 1)
- self.hists.add("Theta_M", 200, 0, 0.02)
+ self.hists.add("Theta_M", 20, 0, 0.001)
self.hists.add("Phi_M", 16, -3.2, 3.2)
- self.hists.add("Mass_M", 200, 0, 2)
+ self.hists.add("Mass_M", 2000, 0, 2)
self.hists.add("ThetaVsP_M", arrayX = tbins, arrayY = pbins)
# Vertex
@@ -163,7 +163,7 @@ class EvgenValidation(EvgenAnalysisAlg):
llp_vtx = None
for i, p in enumerate(evt.particles):
- print("--- ", i)
+ #print("--- ", i)
p.print()
self.fillDaughter(p)
diff --git a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
index 382ad439e5ec1a7d31699c482fd0f17758f18f9c..e361240eb07afde038cfbbda09bf3888521997c4 100644
--- a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
+++ b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
@@ -15,6 +15,8 @@ def getNEvents(fname, maxEvents):
if __name__ == "__main__":
import sys
+
+ doShiftLOS = True
from AthenaCommon.Logging import log
from AthenaCommon.Constants import DEBUG
@@ -49,7 +51,20 @@ if __name__ == "__main__":
cfg = MainServicesCfg(ConfigFlags)
from HEPMCReader.HepMCReaderConfig import HepMCReaderCfg
- cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ if doShiftLOS:
+ cfg.merge(HepMCReaderCfg(ConfigFlags, McEventKey = "BeamTruthEvent_ATLASCoord"))
+ else:
+ cfg.merge(HepMCReaderCfg(ConfigFlags))
+
+ if doShiftLOS:
+ import McParticleEvent.Pythonizations
+ from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg
+ cfg.merge(ShiftLOSCfg(ConfigFlags,
+ xcross = 0,
+ ycross = -150e-6,
+ xshift = 0,
+ yshift = 12))
from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
cfg.merge(McEventSelectorCfg(ConfigFlags))
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index f02f4311f69d07a29a2de4bf26ac2ee93ea13206..3d84bb15354e0bd1ab8bf47eb1a380b27d08939a 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -9,7 +9,7 @@ class ForeseeGenerator(object):
Generate LLP particles within FASER acceptance from FORESEE
"""
- def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, outdir = None, path = '.', randomSeed = 12345):
+ def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, outdir = None, path = '.', randomSeed = 12345, t0 = 0, notime = False, suffix = ""):
self.modelname = modelname
self.energy = energy
@@ -19,8 +19,13 @@ class ForeseeGenerator(object):
self.daughter2_pid = daughter2_pid
self.outdir = outdir
self.path = path
- self.version = 2 # Forsee "version": 2 is in testing
+ self.version = 3 # Forsee "version"
self.seed = randomSeed
+ self.t0 = t0
+ self.notime = notime
+ self.suffix = f"_{suffix}" if suffix else ""
+ self.nbinsample = 1
+ self.nevents = 0
# Set decay mode ...
@@ -42,11 +47,17 @@ class ForeseeGenerator(object):
else:
self.foresee = Foresee(path = self.path)
+ # Generate 6 cm high to account for translation from ATLAS to FASER coord. system
# TODO: relax this a bit as daughters may enter even if mother doesn't
- self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1",
+# 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
+
+ self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.06)**2)< 0.1",
channels=[self.mode], distance=480, length=1.5 ,
luminosity=1/1000.) # 1 pb-1
+
# Set model ...
if self.version == 1:
self.model = Model(self.modelname)
@@ -64,6 +75,8 @@ class ForeseeGenerator(object):
def darkphoton(self):
+ self.nbinsample = 100 # resample bins to help with asymmetric detector
+
# Production modes
self.model.add_production_2bodydecay(
pid0 = "111",
@@ -71,7 +84,7 @@ class ForeseeGenerator(object):
br = "2.*0.99 * coupling**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
generator = "EPOSLHC",
energy = self.energy,
- nsample = 10)
+ nsample = 100)
self.model.add_production_2bodydecay(
pid0 = "221",
@@ -79,7 +92,7 @@ class ForeseeGenerator(object):
br = "2.*0.39 * coupling**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
generator = "EPOSLHC",
energy = self.energy,
- nsample = 10)
+ nsample = 100)
self.model.add_production_mixing(
pid = "113",
@@ -152,13 +165,15 @@ class ForeseeGenerator(object):
def alp_W(self):
+ self.nbinsample = 100 # resample bins to help with smoothness
+
self.model.add_production_2bodydecay(
pid0 = "5",
pid1 = "321",
br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
generator = "Pythia8",
energy = self.energy,
- nsample = 20, # Vary over phi and theta
+ nsample = 500, # Vary over phi and theta -> increase number to improve smoothness for B
)
self.model.add_production_2bodydecay(
@@ -167,7 +182,7 @@ class ForeseeGenerator(object):
br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
generator = "Pythia8",
energy = self.energy,
- nsample = 20,
+ nsample = 500,
)
self.model.add_production_2bodydecay(
@@ -176,7 +191,7 @@ class ForeseeGenerator(object):
br = "4.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
generator = "EPOSLHC",
energy = self.energy,
- nsample = 10,
+ nsample = 50,
)
self.model.add_production_2bodydecay(
@@ -185,9 +200,19 @@ class ForeseeGenerator(object):
br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
generator = "EPOSLHC",
energy = self.energy,
- nsample = 10,
+ nsample = 50,
+ )
+
+ self.model.add_production_2bodydecay(
+ pid0 = "-321",
+ pid1 = "211",
+ br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
+ generator = "EPOSLHC",
+ energy = self.energy,
+ nsample = 50,
)
+
if self.version == 1:
self.model.set_br_1d(
modes = [self.mode],
@@ -229,12 +254,26 @@ class ForeseeGenerator(object):
# Get list of events within detector
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))
+ if self.version >= 3:
+ coups, ctaus, nevents, momenta, weights = output
+ fmomenta = flatten(momenta)
+ fweights = flatten(weights)
+ fenergies = [p.e for p in fmomenta] # Now in MeV ?
+ fthetas = [p.pt/p.pz for p in fmomenta]
- # Return energy (converting to MeV), theta and weights
- return [[e*1000 for e in flatten(energies)], flatten(thetas), flatten(weights)]
+ #self.plot(fthetas, fenergies, fweights)
+
+ self.nevents = nevents
+
+ # Return energy, theta and weights
+ return [fenergies, fthetas, fweights]
+ else:
+ coups, ctaus, nsigs, energies, weights, thetas = output
+ #self.plot(flatten(thetas), flatten(energies), flatten(weights))
+
+ # Return energy (converting to MeV), theta and weights
+ return [[e*1000 for e in flatten(energies)], flatten(thetas), flatten(weights)]
def plot(self, thetas, energies, weights):
# Plot the results in Forsee format
@@ -271,6 +310,13 @@ class ForeseeGenerator(object):
ax.set_ylim(pmin, pmax)
plt.savefig(f"{self.modelname}_m{self.mass}_acc.png")
+
+ def events(self, lumi, f):
+ f.write(f"{self.mass*1000} {self.couplings[0]} {self.nevents[0] * 1000 * lumi}\n")
+ print(f"{self.mass*1000} {self.couplings[0]} {self.nevents[0] * 1000 * lumi}")
+ return
+
+
def write(self):
# Write LLP results to a file
@@ -308,16 +354,19 @@ class ForeseeGenerator(object):
elif not os.path.exists(self.outdir):
os.mkdir(self.outdir)
- filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.hepmc"
+ filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]:.1e}to_{self.daughter1_pid}_{self.daughter2_pid}{self.suffix}.hepmc"
- self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = -1.5, seed = self.seed, decaychannels = [self.mode])
+ _, weights, _ = self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = -1.5, seed = self.seed, decaychannels = [self.mode],
+ notime = self.notime, t0 = self.t0, nsample = self.nbinsample, return_data = True)
cfgname = f"{self.foresee.dirpath}/Models/{self.modelname}/" + filename.replace(".hepmc", ".cfg")
print(f"save config to file: {cfgname}")
with open(cfgname, "w") as f:
f.write(" ".join(sys.argv))
- return
+ filename = f"{self.foresee.dirpath}/Models/{self.modelname}/{filename}"
+
+ return filename, weights
def setup_foresee(path):
@@ -372,22 +421,86 @@ def parse_couplings(data, write_hepMC = False):
return couplings
+def production(args):
+
+ import json
+ import subprocess as proc
+
+ args.hepmc = True
+
+ with open(args.model) as f:
+ data = json.load(f)
+
+ xsect = {}
+ nrun = data["nrun"]
+
+ if args.outdir:
+ outdir = args.outdir + "/" + data["name"]
+ if not os.path.exists(outdir):
+ os.mkdir(outdir)
+ else:
+ outdir = "."
+
+ for coup, masses in data["samples"].items():
+ couplings = parse_couplings([coup], args.hepmc)
+
+ for m in masses:
+ print(f">>> {nrun}: {data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV")
+
+ f = ForeseeGenerator(data["model"], args.Ecom, m, couplings, data["pid1"], data["pid2"], path = args.path, randomSeed = args.randomSeed, t0 = args.t0, notime = args.notime, suffix = data["name"])
+
+ with open("events.txt", "a") as fout:
+ f.events(40, fout)
+ #continue
+
+ if args.hepmc:
+ filename, weights = f.write_hepmc(data["nevents"])
+ else:
+ f.write()
+
+ if args.outdir:
+ p = proc.run(["cp", filename, outdir + "/" + filename.split("/")[-1]])
+ if not p.returncode:
+ proc.run(["rm", filename])
+
+ xsect[nrun] = {
+ "Description" : f"{data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV",
+ "Mass (GeV)" : m,
+ "Coupling" : f"{couplings[0]:e}",
+ "Cross-section (pb)" : sum(weights[0]), # / data["nevents"],
+ "ECom (TeV)" : float(args.Ecom)
+ }
+
+ print(xsect)
+ with open (outdir + "/cross-sections.json", "w") as f:
+ json.dump(xsect, f, indent = 4)
+
+
+ nrun += 1
+
+
+
+ return
+
if __name__ == "__main__":
import argparse, sys
parser = argparse.ArgumentParser(description="Run FORSEE generation")
parser.add_argument("model", help = "Name of foresee model")
- parser.add_argument("--mass", "-m", required = True, type = float, help = "Mass of mother [GeV]")
- parser.add_argument("--couplings", "-c", required = True, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
- parser.add_argument("--pid1", required = True, type = int, help = "PID of daughter 1")
+ parser.add_argument("--mass", "-m", required = False, type = float, help = "Mass of mother [GeV]")
+ parser.add_argument("--couplings", "-c", required = False, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
+ parser.add_argument("--pid1", default = None, required = False, type = int, help = "PID of daughter 1")
parser.add_argument("--pid2", default = None, type = int, help = "PID of daughter 2 (if not set then will be -PID1)")
parser.add_argument("--Ecom", default = "14", help = "Center of mass energy [TeV]")
parser.add_argument("--outdir", "-o", default = None, help = "Output path")
parser.add_argument("--path", default = ".", help = "Path to foresee installation")
parser.add_argument("--hepmc", action = "store_true", help = "Write HepMC events")
parser.add_argument("--nevents", "-n", default = 10, type = int, help = "Number of HepMC events ")
- parser.add_argument("--randomSeed", "-s", default = 1234, type = int, help = "Random seed for HepMC generation")
+ parser.add_argument("--randomSeed", "-s", default = 1234, type = int, help = "Random seed for HepMC generation")
+ parser.add_argument("--t0", "-t", default = 0, type = int, help = "Time offset for start of decay volume")
+ parser.add_argument("--notime", action = "store_true", help = "Set all vertex times to 0 rather than calculating from start of decay volume")
+ parser.add_argument("--suffix", default = "", help = "Filename suffix")
args = parser.parse_args()
add_to_python_path(f"{args.path}/src")
@@ -395,22 +508,26 @@ if __name__ == "__main__":
from foresee import Foresee, Model, Utility
# Create PIDs
- if args.pid2 is None:
+ if args.pid2 is None and args.pid1 is not None:
args.pid2 = -args.pid1
-
- couplings = parse_couplings(args.couplings, args.hepmc)
-
- print(f"Generating {args.model} events at Ecom = {args.Ecom}")
- print(f" mother mass = {args.mass} GeV")
- print(f" decay = {args.pid1} {args.pid2}")
- print(f" couplings = {couplings}")
- f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, outdir = args.outdir, path = args.path, randomSeed = args.randomSeed)
+ if args.model.endswith(".json"):
+ production(args)
+ else:
+
+ couplings = parse_couplings(args.couplings, args.hepmc)
- if args.hepmc:
- f.write_hepmc(args.nevents)
- else:
- f.write()
+ print(f"Generating {args.model} events at Ecom = {args.Ecom}")
+ print(f" mother mass = {args.mass} GeV")
+ print(f" decay = {args.pid1} {args.pid2}")
+ print(f" couplings = {couplings}")
+
+ f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, outdir = args.outdir, path = args.path, randomSeed = args.randomSeed, t0 = args.t0, notime = args.notime, suffix = args.suffix)
+
+ if args.hepmc:
+ f.write_hepmc(args.nevents)
+ else:
+ f.write()
diff --git a/Generators/ForeseeGenerator/share/plot_validation.py b/Generators/ForeseeGenerator/share/plot_validation.py
index 8429978bc6a091abd853456d1a10fa3d33ab200d..d4ff9c8b2cb9274a7aade2e8b0a0f6ceb6bf21e9 100644
--- a/Generators/ForeseeGenerator/share/plot_validation.py
+++ b/Generators/ForeseeGenerator/share/plot_validation.py
@@ -99,7 +99,7 @@ if __name__ == "__main__":
plotn(f, args, config, 3, 2, f"daug{i}")
config = [Hist("P_M", logy = True, xtitle = "p^{0} [GeV]", ndiv = 5, xlo = 0, xhi = 10000, r=10),
- Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4, r = 10),
+ Hist("Theta_M", xtitle = "#theta [rad]", ndiv = -4, r = 1),
Hist("Mass_M", xtitle = "m^{0} [GeV]", xlo = 0, xhi = 1., ndiv = 5),
Hist("Pt_M", logy = True, xtitle = "p_{T}^{0} [GeV]", ndiv = 10, r = 50),
Hist("Phi_M", xtitle = "#phi [rad]", r = 2),
@@ -111,31 +111,31 @@ if __name__ == "__main__":
plotn(f, args, Hist("PIDs", xtitle="PDG Id"), 1, 1, "pid")
-# config = [Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
-# Hist("ThetaVsP_d0", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
-# Hist("ThetaVsP_d1", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz")
-# ]
-#
-# plotn(f, args, config, 2, 2, "twod")
-
- config = [Hist("Vtx_X_LLP", xtitle = "x [mm]", r = 5),
- Hist("Vtx_Y_LLP", xtitle = "y [mm]", r = 5),
- Hist("Vtx_Z_LLP", xtitle = "z [mm]", r = 5, ndiv = 5),
- Hist("Vtx_XY_LLP", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
- Hist("Vtx_R_LLP", xtitle = "r [mm]", r = 5, ndiv = 5)
+ config = [Hist("ThetaVsP_M", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
+ Hist("ThetaVsP_d0", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz"),
+ Hist("ThetaVsP_d1", ytitle = "p^{0} [GeV]", xtitle = "#theta [rad]", logx = True, logy = True, d = "colz")
+ ]
+
+ plotn(f, args, config, 2, 2, "twod")
+
+ config = [Hist("Vtx_X_LLP", xtitle = "x [mm]", r = 1),
+ Hist("Vtx_Y_LLP", xtitle = "y [mm]", r = 1),
+ Hist("Vtx_Z_LLP", xtitle = "z [mm]", r = 1, ndiv = 5),
+ Hist("Vtx_XY_LLP", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (1,1)),
+ Hist("Vtx_R_LLP", xtitle = "r [mm]", r = 1, ndiv = 5)
]
plotn(f, args, config, 3, 2, "vtx_llp")
- config = [Hist("Vtx_X_All", xtitle = "x [mm]", r = 5),
- Hist("Vtx_Y_All", xtitle = "y [mm]", r = 5),
- Hist("Vtx_Z_All", xtitle = "z [mm]", r = 5, ndiv = 5),
- Hist("Vtx_XY_All", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
- Hist("Vtx_R_All", xtitle = "r [mm]", r = 5, ndiv = 5)
- ]
-
- plotn(f, args, config, 3, 2, "vtx_all")
+# config = [Hist("Vtx_X_All", xtitle = "x [mm]", r = 5),
+# Hist("Vtx_Y_All", xtitle = "y [mm]", r = 5),
+# Hist("Vtx_Z_All", xtitle = "z [mm]", r = 5, ndiv = 5),
+# Hist("Vtx_XY_All", xtitle = "x [mm]", ytitle = "y [mm]", d = "colz", r = (5,5)),
+# Hist("Vtx_R_All", xtitle = "r [mm]", r = 5, ndiv = 5)
+# ]
+#
+# plotn(f, args, config, 3, 2, "vtx_all")
# config = [Hist("Vtx_X", xtitle = "x [mm]", r = 5),
diff --git a/Generators/ForeseeGenerator/share/validate_grid.py b/Generators/ForeseeGenerator/share/validate_grid.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3d5302dff75a44f6cc2a77497f625e6659706db
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/validate_grid.py
@@ -0,0 +1,36 @@
+import json
+import subprocess as proc
+import os
+
+grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
+
+with open(grid) as f:
+ data = json.load(f)
+
+name = "DarkPhotonGrid2"
+path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
+ecom = 14
+pid1 = 11
+pid2 = -11
+
+for coup, masses in data["samples"].items():
+ for m in masses:
+ infile = f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):e}to_{pid1}_{pid2}_{name}.hepmc"
+ outfile = infile.replace(".hepmc", ".EVNT.pool.root")
+ valfile = infile.replace(".hepmc", ".VAL.pool.root")
+ valdir = infile.replace(".hepmc", "")
+
+ #if str(coup) in coup_fix:
+ # infile = infile.replace(str(coup), coup_fix[str(coup)])
+
+ cmd = f"python ../calypso/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py --filesInput {infile} Output.EVNTFileName={outfile}"
+ print(">>>", cmd, "\n")
+ os.system(cmd)
+
+ cmd2 = f"python ../calypso/Generators/ForeseeGenerator/python/Validate.py -d 2 -m -n 10000 -o {valfile} {outfile}"
+ print(">>>", cmd2, "\n")
+ os.system(cmd2)
+
+ cmd3 = f"python ../calypso/Generators/ForeseeGenerator/share/plot_validation.py -d 2 -o {valdir} {valfile}"
+ print(">>>", cmd3, "\n")
+ os.system(cmd3)