Skip to content
Snippets Groups Projects
Commit 26b84127 authored by Carl Gwilliam's avatar Carl Gwilliam
Browse files

Update generator scripts for B-L and add grids for ALPS models

parent 9c8341c8
No related branches found
No related tags found
2 merge requests!354Update generator scripts for B-L and add grids for ALPS models,!351Merging Muon code to create new Database tag
Showing
with 232 additions and 10 deletions
......@@ -109,13 +109,13 @@ class EvgenValidation(EvgenAnalysisAlg):
# Vertex
self.hists.add("Vtx_X_LLP", 50, -100, 100)
self.hists.add("Vtx_Y_LLP", 50, -100, 100)
self.hists.add("Vtx_Z_LLP", 500, -5000, 0)
self.hists.add("Vtx_Z_LLP", 500, -5000, 5000)
self.hists.add("Vtx_R_LLP", 20, 0, 200)
self.hists.add("Vtx_XY_LLP", 50, -100, 100, 50, -100, 100)
self.hists.add("Vtx_X_All", 50, -100, 100)
self.hists.add("Vtx_Y_All", 50, -100, 100)
self.hists.add("Vtx_Z_All", 500, -5000, 0)
self.hists.add("Vtx_Z_All", 500, -5000, 5000)
self.hists.add("Vtx_R_All", 20, 0, 200)
self.hists.add("Vtx_XY_All", 50, -100, 100, 50, -100, 100)
......
import os
import os, sys
import numpy as np
import matplotlib.pyplot as plt
......@@ -53,8 +53,17 @@ class ForeseeGenerator(object):
# 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.065)**2)< 0.1",
channels=[self.mode], distance=480, length=1.5 ,
# Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
# Relax r to 11 cm to acount for any mismatches with calypso and for duaghter entering if mother doesn't
if "ALP" in self.modelname:
# Since not relying on tracks then extend length to include spectrometer (magents + tracking stations)
length = 1.5 + 2 + 0.5 # m
else:
length = 1.5 # m
self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.11",
channels=[self.mode], distance=480, length=length ,
luminosity=1/1000.) # 1 pb-1
......@@ -68,6 +77,10 @@ class ForeseeGenerator(object):
self.data = self.darkphoton()
elif self.modelname == "ALP-W":
self.data = self.alp_W()
elif self.modelname == "ALP-photon":
self.data = self.alp_photon()
elif self.modelname == "U(1)B-L":
self.data = self.bminusl()
else:
sys.exit(f"Unknown model {self.modelname}")
......@@ -227,6 +240,84 @@ class ForeseeGenerator(object):
)
return self.decays()
def alp_photon(self):
if self.version == 1:
sys.exit("ALP-photon model not available in this version of FORESEE")
self.nbinsample = 100 # resample bins to help with smoothness
masses = [float(x) for x in ['{:0.4e}'.format(m) for m in np.logspace(-2,0,20+1)]]
self.model.add_production_direct(
label = "Prim",
energy = self.energy,
coupling_ref = 1,
masses=masses,
)
self.model.set_br_1d(
modes = ["gamma_gamma"],
finalstates = [[22, 22]],
filenames=["model/br/gamma_gamma.txt"]
)
return self.decays()
def bminusl(self):
if self.version == 1:
sys.exit("U(1)B-L model not available in this version of FORESEE")
self.nbinsample = 100 # TODO
self.model.add_production_2bodydecay(
pid0 = "111",
pid1 = "22",
br = "2.*0.99 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
generator = "EPOSLHC",
energy = self.energy,
nsample = 100)
self.model.add_production_2bodydecay(
pid0 = "221",
pid1 = "22",
br = "2.*0.25*0.39 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
generator = "EPOSLHC",
energy = self.energy,
nsample = 100)
masses_brem = [
0.01 , 0.0126, 0.0158, 0.02 , 0.0251, 0.0316, 0.0398,
0.0501, 0.0631, 0.0794, 0.1 , 0.1122, 0.1259, 0.1413,
0.1585, 0.1778, 0.1995, 0.2239, 0.2512, 0.2818, 0.3162,
0.3548, 0.3981, 0.4467, 0.5012, 0.5623, 0.6026, 0.631 ,
0.6457, 0.6607, 0.6761, 0.6918, 0.7079, 0.7244, 0.7413,
0.7586, 0.7762, 0.7943, 0.8128, 0.8318, 0.8511, 0.871 ,
0.8913, 0.912 , 0.9333, 0.955 , 0.9772, 1. , 1.122 ,
1.2589, 1.4125, 1.5849, 1.7783, 1.9953, 2.2387, 2.5119,
2.8184, 3.1623, 3.9811, 5.0119, 6.3096, 7.9433, 10.
]
self.model.add_production_direct(
label = "Brem",
energy = self.energy,
condition = "1./0.3**2 * (p.pt<1)",
coupling_ref=1,
masses = masses_brem,
)
self.model.set_br_1d(
modes = ["e_e", "mu_mu", "nu_nu", "Hadrons"],
finalstates = [[11, -11], [13, -13], [12, -12], None],
filenames=["model/br/e_e.txt", "model/br/mu_mu.txt", "model/br/nu_nu.txt", "model/br/Hadrons.txt"],
#filenames=[f"model/br/all.txt"]
)
return self.decays()
def decays(self):
......@@ -468,7 +559,8 @@ def production(args):
"Mass (GeV)" : m,
"Coupling" : f"{couplings[0]:e}",
"Cross-section (pb)" : sum(weights[0]), # / data["nevents"],
"ECom (TeV)" : float(args.Ecom)
"ECom (TeV)" : float(args.Ecom),
"NEvents" : float(data["nevents"])
}
print(xsect)
......
{
"name" : "ALPPhotontGrid13p6_65mm",
"model" : "ALP-photon",
"pid1" : 22,
"pid2" : 22,
"nevents" : 50000,
"nrun" : 110100,
"samples" : {
"7e-05" : [0.06],
"1.5e-04" : [0.06],
"3e-04" : [0.06]
}
}
{
"name" : "DarkPhotonGrid13p6_65mm",
"model" : "DarkPhoton",
"pid1" : 11,
"pid2" : -11,
"nevents" : 50000,
"nrun" : 110021,
"samples" : {
"1e-06" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794],
"5e-06" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
"1e-05" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
"3e-05" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
"5e-05" : [0.01, 0.0251, 0.0501, 0.0631],
"1e-04" : [0.01, 0.0251, 0.0501],
"3e-04" : [0.01, 0.0251]
}
}
{
"name" : "ALPWGrid13p6_65mm_110mm_3p5m_extra",
"model" : "ALP-W",
"pid1" : 22,
"pid2" : 22,
"nevents" : 50000,
"nrun" : 110206,
"samples" : {
"4e-04" : [0.1],
"1e-04" : [0.2],
"5e-05" : [0.1],
"2e-05" : [0.35]
}
}
{
"name" : "ALPPhotonGrid13p6_65mm_110mm_3p5m",
"model" : "ALP-photon",
"pid1" : 22,
"pid2" : 22,
"nevents" : 50000,
"nrun" : 110103,
"samples" : {
"7e-05" : [0.06],
"1.5e-04" : [0.06],
"3e-04" : [0.06]
}
}
{
"name" : "U1B-LGrid13p6_65mm_110mm",
"model" : "U(1)B-L",
"pid1" : 11,
"pid2" : -11,
"nevents" : 50000,
"nrun" : 110300,
"samples" : {
"4e-5" : [0.01],
"3e-5" : [0.01],
"2.4e-5" : [0.014],
"2e-5" : [0.014],
"1.4e-5" : [0.02],
"1.1e-5" : [0.02, 0.024],
"1e-5" : [0.024],
"8e-6" : [0.03],
"7e-6" : [0.03],
"6e-6" : [0.03],
"5e-6" : [0.03, 0.034],
"4e-6" : [0.01, 0.014, 0.02, 0.03, 0.034, 0.04],
"3e-6 " : [0.01, 0.014, 0.02, 0.024, 0.03, 0.034, 0.04],
"2.4e-6" : [ 0.024, 0.03]
}
}
{
"name" : "DarkPhotonGrid13p6_65mm_extra",
"model" : "DarkPhoton",
"pid1" : 11,
"pid2" : -11,
"nevents" : 50000,
"nrun" : 110056,
"samples" : {
"3e-04" : [0.0126, 0.0158],
"1e-04" : [0.0126, 0.0158],
"5e-05" : [ 0.0158, 0.0316, 0.0358],
"3e-05" : [ 0.0316, 0.0358],
"2e-05" : [ 0.0501, 0.0631, 0.0794],
"1e-05" : [0.0126, 0.0158, 0.0316, 0.0358],
"5e-06" : [0.0126, 0.0158, 0.0316, 0.0358]
"3e-06" : [0.01, 0.0126, 0.0158, 0.0251, 0.0316, 0.0358, 0.0501, 0.0631, 0.0794],
}
}
{
"name" : "DarkPhotonGrid13p6_65mm_extra2",
"model" : "DarkPhoton",
"pid1" : 11,
"pid2" : -11,
"nevents" : 50000,
"nrun" : 110090,
"samples" : {
"1.5e-04" : [0.01, 0.0126, 0.0158, 0.0251]
}
}
{
"name" : "DarkPhotonGrid13p6b",
"model" : "DarkPhoton",
"pid1" : 11,
"pid2" : -11,
"nevents" : 50000,
"nrun" : 110056,
"samples" : {
"3e-06" : [0.0398]
}
}
......@@ -3,19 +3,20 @@ import subprocess as proc
from glob import glob
import os, sys
grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
grid = "../calypso/Generators/ForeseeGenerator/share/points_ALP_W.json"
with open(grid) as f:
data = json.load(f)
name = "DarkPhotonGrid13p6_65mm"
name = data["name"]
path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
ecom = 13.6
pid1 = 11
pid2 = -11
pid1 = data["pid1"]
pid2 = data["pid2"]
for coup, masses in data["samples"].items():
for m in masses:
print(coup, m, f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}*.hepmc")
files = glob(f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}*.hepmc")
if len(files) != 1:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment