Skip to content
Snippets Groups Projects
Commit 422fea3e authored by FaserMC's avatar FaserMC
Browse files

Add Carl's DIF update

parents ae897693 5c1f38f5
Branches
Tags
1 merge request!333Physics Ntuple update
......@@ -47,6 +47,12 @@ def faser_pgparser():
parser.add_argument("--zpos", default=None, type=float,
help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
parser.add_argument("--pidd1", default=None, type=int,
help="Specify PDG ID of daugther 1 for DIF generator")
parser.add_argument("--pidd2", default=None, type=int,
help="Specify PDG ID of daugther 2 for DIF generator")
parser.add_argument("--sampler", default="log",
help="Specify energy sampling (log, lin, const, hist, hist2D)")
parser.add_argument("--hist_name", default="log",
......
......@@ -101,8 +101,11 @@ if __name__ == '__main__':
from AthenaCommon.PhysicalConstants import pi
if isinstance(args.pid, list):
# Note args.pid is a list, must make this a set for ParticleGun
pidarg = set(args.pid)
# Note args.pid is a list, must make this a set for ParticleGun if have > 1
if len(args.pid) > 1:
pidarg = set(args.pid)
else:
pidarg = args.pid[0]
else:
# Just pass a single value
pidarg = args.pid
......@@ -112,13 +115,32 @@ if __name__ == '__main__':
# Create the simgun dictionary
# Negative radius gives uniform sampling
# Positive radius gives Gaussian sampling
sg_dict = {
"Generator" : "SingleParticle",
"pid" : pidarg, "mass" : args.mass,
"theta" : PG.GaussianSampler(0, args.angle, oneside = True),
"phi" : [0, 2*pi], "radius" : args.radius,
"randomSeed" : args.outfile
}
if args.pidd1:
if args.pidd2 is None:
args.pidd2 = -args.pidd1
sg_dict = {
"Generator" : "DecayInFlight",
"mother_pid" : pidarg,
"daughter1_pid" : args.pidd1,
"daughter2_pid" : args.pidd2,
"mass" : args.mass,
"theta" : PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
"phi" : [0, 2*pi], "radius" : args.radius,
"randomSeed" : args.outfile
}
else:
sg_dict = {
"Generator" : "SingleParticle",
"pid" : pidarg, "mass" : args.mass,
"theta" : PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
"phi" : [0, 2*pi], "radius" : args.radius,
"randomSeed" : args.outfile
}
# -1000 is safely upstream of detector (to be checked)
# Note zpos is in mm!
......@@ -132,6 +154,12 @@ if __name__ == '__main__':
sg_dict["energy"] = PG.LogSampler(args.minE*GeV, args.maxE*GeV)
elif args.sampler == "const":
sg_dict["energy"] = PG.ConstSampler(args.maxE*GeV)
elif args.sampler == "hist":
fname, hname = args.hist_name.split(":")
sg_dict["energy"] = PG.TH1Sampler(fname, hname)
elif args.sampler == "hist2D":
fname, hname = args.hist_name.split(":")
sg_dict["energy"] = PG.TH2Sampler(fname, hname)
else:
print(f"Sampler {args.sampler} not known!")
sys.exit(1)
......
......@@ -93,9 +93,8 @@ class DIFSampler(PG.ParticleSampler):
# def __init__(self, daughter1_pid = 11, daughter2_pid = -11, mother_pid = None, mother_mom = PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (10*MeV)**2),0,10*MeV,0),mother_pos = CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0)):
def __init__(self, daughter1_pid=13, daughter2_pid=-13, mother_pid=None,
mother_mom=PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (500*MeV)**2), [0, 0.0002], 500*MeV, [0, 2*pi]),
my_z_position=-1500):
# mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=CylinderSampler([0, 100**2], [0, 2*pi], my_z_position, 0))
mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=mother_pos)
self.daughter1 = self.particle(daughter1_pid)
self.daughter2 = self.particle(daughter2_pid)
......
......@@ -9,7 +9,7 @@ import ParticleGun as PG
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
# from AthenaCommon.Constants import VERBOSE, INFO
from AthenaCommon.SystemOfUnits import TeV
from AthenaCommon.SystemOfUnits import TeV, GeV, MeV
from AthenaCommon.PhysicalConstants import pi
### add radial pos sampler ### with gaussian beam implemented
......@@ -144,7 +144,31 @@ def FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs) :
pg = cfg.getPrimary()
from DIFGenerator import DIFSampler
pg.sampler = DIFSampler(**kwargs)
if "radius" in kwargs:
kwargs["mother_pos"] = RadialPosSampler(x = kwargs.setdefault("x", 0.0),
y = kwargs.setdefault("y", 0.0),
z = kwargs.setdefault("z", -3750.0),
r = kwargs.setdefault("radius", 1.0),
t = kwargs.setdefault("t", 0.0) )
else:
from DIFGenerator.DIFSampler import CylinderSampler
kwargs["mother_pos"] = CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)
if not "mother_mom" in kwargs:
kwargs["mother_mom"] = PG.EThetaMPhiSampler(kwargs.setdefault("energy", ((1*TeV)**2 + (500*MeV)**2)**0.5),
kwargs.setdefault("theta", [0, 0.0002]),
kwargs.setdefault("mass", 500*MeV),
kwargs.setdefault("phi", [0, 2*pi]) )
pg.sampler = DIFSampler(kwargs.setdefault("daughter1_pid", 13),
kwargs.setdefault("daughter2_pid", -13),
kwargs.setdefault("mother_pid", None),
kwargs["mother_mom"],
kwargs["mother_pos"] )
return cfg
......@@ -190,9 +214,12 @@ def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
return cfg
def FaserParticleGunCfg(ConfigFlags) :
# generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
# generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
kwargs = ConfigFlags.Sim.Gun
del kwargs["Generator"]
if generator == "SingleEcalParticle" :
return FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs)
elif generator == "Cosmics" :
......
......@@ -47,13 +47,13 @@ class ForeseeGenerator(object):
else:
self.foresee = Foresee(path = self.path)
# Generate 6 cm high to account for translation from ATLAS to FASER coord. system
# Generate 6.5 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",
# 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",
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 ,
luminosity=1/1000.) # 1 pb-1
......
import json
import subprocess as proc
import os
from glob import glob
import os, sys
grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
with open(grid) as f:
data = json.load(f)
name = "DarkPhotonGrid13p6"
name = "DarkPhotonGrid13p6_65mm"
path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
ecom = 13.6
pid1 = 11
......@@ -15,7 +16,12 @@ 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):.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:
continue
infile = files[0]
outfile = infile.replace(".hepmc", ".EVNT.pool.root")
valfile = infile.replace(".hepmc", ".VAL.pool.root")
valdir = infile.replace(".hepmc", "")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment