diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
index 2854fd9932e0b6e561020984bac57b8b1f39b86c..a99c51eecdc4a813901085de67d792be40f73cb6 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
@@ -131,7 +131,7 @@ if len(args.short) > 0:
# more files later. --complete overrides this if wanted
# Build output filename
-if seglo == 0 and (seghi+1) == len(dirlist) and args.complete: # Full run
+if args.complete and (seglo == 0) and ((seghi+1) == len(dirlist)): # Full run
outfile = f"FaserMC-{short}-{run}"
elif seglo == seghi: # Single segment
outfile = f"FaserMC-{short}-{run}-{seglo:05}"
diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index ff8f1d3792a666bf92dca51b16bf7adf58c7500d..c98d348aea5ccc2e2b3c1459ef793bd41bb8c3e9 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -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",
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index cc1628efdb9245e32b7abbbc9ca60a60b2151fad..d23e879d29feefe759d8dcb5a617cfea3bebb1c2 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -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)
diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
index 5f1ee06615cd822c911cc67d0acc7659662e78dc..f75b144a1bf2d0a725692dbaa82be34ad788cd5c 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
+++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
@@ -11,7 +11,7 @@
#
# Monte Carlo options:
# --isMC - needed for MC reco
-# --digiTag <tag> - override MC reco tag for calo gain (matches digi tag)
+# --caloTag <tag> - override MC reco tag for calo gain (to match digi tag)
#
# file_path - full file name (with path)
# release_directory - optional path to release install directory (default pwd)
@@ -55,7 +55,7 @@ do
ismc=1
shift;;
- --digiTag)
+ --caloTag)
echo "Override calo digi tag with $2"
gainstr="--MC_calibTag $2"
shift;
diff --git a/Generators/DIFGenerator/python/DIFSampler.py b/Generators/DIFGenerator/python/DIFSampler.py
index 8ec22459aba93a430093c4bd7b2a3f13e2d0f206..007830b63c82ec06f19af6fbd3169cfc10ce8e2e 100644
--- a/Generators/DIFGenerator/python/DIFSampler.py
+++ b/Generators/DIFGenerator/python/DIFSampler.py
@@ -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)
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index ac30494a4815f5bb06a5d67408b4c0551297a63b..d2e82a88d3fd6229b2e3ad33932ead7d5d650720 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -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" :
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 0c3096f8abe6a640391999613442ce1a267f4aca..6ef5d1d0377a147246b97b9e203fcc224355bae2 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -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
diff --git a/Generators/ForeseeGenerator/share/validate_grid.py b/Generators/ForeseeGenerator/share/validate_grid.py
index cf01468483c71f68b3687ec2a7a17a031fe84d7a..88edb1194ee6ecebe37464906b26a03afde7b1da 100644
--- a/Generators/ForeseeGenerator/share/validate_grid.py
+++ b/Generators/ForeseeGenerator/share/validate_grid.py
@@ -1,13 +1,14 @@
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", "")
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index 5ae87306e068ac7ebfde95f0c29aed5dba283f43..a22d1edbb6802f298044c63e3381d599ce35c978 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -60,7 +60,7 @@ filelist = []
# If this is a directory, need to create file list
if filepath.is_dir():
- # Use expected data pattern to find files
+ # Use expected data pattern to find files (data only)
runstr = filepath.stem
# Make list of segments to search for
@@ -92,7 +92,7 @@ if filepath.is_dir():
for seg in seglist:
if args.isMC:
- searchstr = f"FaserMC-*-{runstr}-{seg}-*xAOD.root"
+ searchstr = f"FaserMC-*-{seg}-*xAOD.root"
else:
searchstr = f"Faser-Physics-{runstr}-{seg}-*xAOD.root"
@@ -114,10 +114,12 @@ if filepath.is_dir():
filelist.append(filestr)
# End of loop over segments
+ # Parse name to create outfile
firstfile = Path(filelist[0])
firststem = str(firstfile.stem)
firstfaser = firststem.split('-')[0]
firstshort = firststem.split('-')[1]
+ runstr = firststem.split('-')[2]
firstseg = firststem.split('-')[3]
if args.merge > 1:
firstseg2 = firststem.split('-')[4]
@@ -128,13 +130,21 @@ if filepath.is_dir():
if args.merge > 1:
lastseg = laststem.split('-')[4]
+ print(f"Faser = {firstfaser}")
+ print(f"Short = {firstshort}")
+ print(f"Run = {runstr}")
+ print(f"First = {firstseg}")
+ print(f"Last = {lastseg}")
+ print(f"Args = {args.tag}")
+
# Find any tags
tagstr = firststem.replace(f"{firstfaser}-{firstshort}-{runstr}-{firstseg}", "")
if args.merge > 1:
tagstr = tagstr.replace(f"-{firstseg2}", "")
tagstr = tagstr.replace("-xAOD", "")
- print(f"Tag = {tagstr}")
+
+ print(f"Tag = {tagstr}")
# Build output name
outfile = f"{firstfaser}-{firstshort}-{runstr}-{firstseg}-{lastseg}"
@@ -148,6 +158,7 @@ if filepath.is_dir():
outfile += "-PHYS.root"
+
# If this is a single file, just process that
# Could be a url, so don't check if this is a file
else:
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 744606b00d5cfef555bfc4992589ff40943f9fc7..e51d97a6a36f9a1b13c9e343200911dc9506cbf2 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -612,6 +612,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
// load in calibrated calo container
+ if (m_doCalib) {
SG::ReadHandle<xAOD::CalorimeterHitContainer> ecalCalibratedContainer { m_ecalCalibratedContainer, ctx };
ATH_CHECK(ecalCalibratedContainer.isValid());
for (auto hit : *ecalCalibratedContainer) {
@@ -651,6 +652,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("Calibrated preshower: ch is " << ch << ", edep is " << hit->E_dep() << ", E_EM is " << hit->E_EM() << ", Nmip is " << hit->Nmip() << ", fit_to_raw_ratio is " << hit->fit_to_raw_ratio());
}
+ }
// process all waveeform data fro all scintillator and calorimeter channels
SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
@@ -776,8 +778,16 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// loop over all reconstructed tracks and use only the tracks that have hits in all three tracking stations (excludes IFT)
// store track parameters at most upstream measurement and at most downstream measurement
// extrapolate track to all scintillator positions and store extrapolated position and angle
- SG::ReadHandle<TrackCollection> trackCollection {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
+ SG::ReadHandle<TrackCollection> trackCollection;
+ if (m_useIFT) {
+ SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that excludes IFT
+ trackCollection = tc;
+ } else {
+ SG::ReadHandle<TrackCollection> tc {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
+ trackCollection = tc;
+ }
ATH_CHECK(trackCollection.isValid());
+
for (const Trk::Track* track : *trackCollection)
{
if (track == nullptr) continue;
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 4e694c21defbc4e284aecd5bb6e40cf4e2360be7..4d1ec93b8f073a68cb31cdece4456e36cb92eee0 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -101,6 +101,8 @@ private:
const PreshowerID* m_preshowerHelper;
const EcalID* m_ecalHelper;
+ BooleanProperty m_useIFT { this, "UseIFT", false, "Use IFT tracks" };
+ BooleanProperty m_doCalib { this, "DoCalib", true, "Fill calibrated calorimeter quantities" };
BooleanProperty m_doBlinding { this, "DoBlinding", true, "Blinding will not output events with Calo signal > 10 GeV e-" };
BooleanProperty m_useFlukaWeights { this, "UseFlukaWeights", false, "Flag to weight events according to value stored in HepMC::GenEvent" };
BooleanProperty m_useGenieWeights { this, "UseGenieWeights", false, "Flag to weight events according to Genie luminosity" };