diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 9f62fd57b12fa450f6e2218b841501ca2351cf2a..c3f5313c0e52ddc68009ac2d6d3e41103b103b8a 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -6,6 +6,7 @@
#include <map>
#include <utility>
+#include <cmath>
CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name,
ISvcLocator* pSvcLocator)
@@ -99,19 +100,31 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
// Create structure to store pulse for each channel
std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
+ // Sum energy for each channel (i.e. identifier)
+ std::map<unsigned int, float> esum;
+ for (const auto& hit : *caloHitHandle) {
+ esum[hit.identify()] += hit.energyLoss();
+ }
+
+ // Loop over time samples
for (const auto& tk : m_timekernel) {
std::map<unsigned int, float> counts;
// Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
- for (const auto& hit : *caloHitHandle) {
- counts[hit.identify()] += tk * hit.energyLoss();
+ //for (const auto& hit : *caloHitHandle) {
+ // counts[hit.identify()] += tk * hit.energyLoss();
+ //}
+
+ // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+ for (const auto& e : esum) {
+ counts[e.first] = tk * e.second;
}
// Subtract count from basleine and add result to correct waveform vector
for (const auto& c : counts) {
- unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = baseline - c.second;
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = std::round(baseline - c.second);
if (value < 0) {
ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi.py b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
index 0e48e1e27dee3b1e17ecddd3fcdef4c17fc0936e..0da1361cec9191f96014543a68aaca777e0acff8 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
@@ -25,6 +25,8 @@ parser.add_argument("-g", "--geom", default="TI12MC",
help="Specify geometry (default: TI12MC, alt: TestBeamMC)")
parser.add_argument("-t", "--tag", default="",
help="Specify digi tag (to append to output filename)")
+parser.add_argument("--subtractTime", type=float, default=-999.,
+ help="Subtract time parameter for SCT RDOs")
parser.add_argument("--digiTag", default="",
help="Specify tag for waveform digi folder")
parser.add_argument("--short", default="",
@@ -161,6 +163,10 @@ acc.merge(FaserGeometryCfg(ConfigFlags))
from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_DigitizationCfg
acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
+# Set the time offset for SCT RDOs
+pualg = acc.getEventAlgo("StandardPileUpToolsAlg")
+pualg.PileUpTools['FaserSCT_DigitizationTool'].SurfaceChargesGenerator.SubtractTime = args.subtractTime
+
# Pass something to set folder tag
from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
index ec595160f8766d3737aef15c53082c970c7e9062..b945ffbbd19c024d549d0d0f84fdfca768759f68 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
@@ -29,8 +29,12 @@ parser.add_argument("-s", "--slice", type=int, default=0,
help="Specify file slice to produce")
parser.add_argument("-f", "--files", type=int, default=5,
help="Specify number of input files to run in one batch")
+parser.add_argument("--complete", action="store_true",
+ help="This is the complete run, remove segment numbers")
parser.add_argument("-t", "--tag", default="",
help="Specify digi tag (to append to output filename)")
+parser.add_argument("--subtractTime", type=float, default=-999.,
+ help="Subtract time parameter for SCT RDOs")
parser.add_argument("--digiTag", default="",
help="Specify tag for waveform digi folder")
parser.add_argument("--short", default="",
@@ -67,24 +71,35 @@ if len(dirlist) == 0:
print("HITS files available:")
[print(file) for file in dirlist]
+# Get run number for pattern matching
+stem = dirlist[0].stem
+spl = stem.split('-')
+short = spl[1]
+run = spl[2]
+
for seg in seglist:
- # Assume these are in numerical order from 0
- if seg >= len(dirlist):
- print(f"Requested file segment {seg} but only {len(dirlist)} files found")
- if args.partial:
- break
- else:
- sys.exit(1) # Abort this job
+ # No longer assume all HITS files are complete
+ #if seg >= len(dirlist):
+ # print(f"Requested file segment {seg} but only {len(dirlist)} files found")
+ # if args.partial:
+ # break
+ # else:
+ # sys.exit(1) # Abort this job
# Check if segment number exists in hits file (this is not perfect)
- segstr = f"{seg:05d}"
- if segstr not in dirlist[seg].name:
- print(f"Segment {segstr} not in file {dirlist[seg].name}!")
+ segstr = f"FaserMC-{short}-{run}-{seg:05d}*HITS.root"
+ matchlist = list(dirpath.glob(segstr))
+ if len(matchlist) == 0:
+ print(f"Segment {segstr} not found!")
if not args.partial: sys.exit(1) # abort
+ continue
+ elif len(matchlist) > 1:
+ print(f"Segment {segstr} matched {len(matchlist)} files!")
+ sys.exit(1) # Should never happen!
else:
- print(f"Segment {segstr} found in file {dirlist[seg]}")
+ print(f"Segment {segstr} found!")
- filelist.append(dirlist[seg])
+ filelist.append(matchlist[0])
if len(filelist) == 0:
# Asked for range that doesn't exist
@@ -114,8 +129,11 @@ seghi = int(spl[3])
if len(args.short) > 0:
short += args.short
+# Shouldn't normally remove segments with MC as we may always add
+# more files later. --complete overrides this if wanted
+
# Build output filename
-if seglo == 0 and (seghi+1) == len(dirlist): # 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}"
@@ -215,6 +233,10 @@ acc.merge(FaserGeometryCfg(ConfigFlags))
from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_DigitizationCfg
acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
+# Set the time offset for SCT RDOs
+pualg = acc.getEventAlgo("StandardPileUpToolsAlg")
+pualg.PileUpTools['FaserSCT_DigitizationTool'].SurfaceChargesGenerator.SubtractTime = args.subtractTime
+
# Pass something to set folder tag
from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
index 7ab8eb3734b2cba6676ecd8651967233565ae6dc..fcf5e2045b192ef913a9fc1ea0b75c8fae4074bd 100644
--- a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi.sh
@@ -2,10 +2,10 @@
# Used with a condor file to submit to vanilla universe
#
# Usage:
-# submit_faser_digi.sh [--highGain] filepath [release_directory] [working_directory]
+# submit_faser_digi.sh filepath [release_directory] [working_directory]
#
# Options:
-# --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
+# --digiTag <tag> - override digitization tag for calo gain
# --geom - geometry setting
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
@@ -34,6 +34,12 @@ do
shift;
shift;;
+ --subtractTime)
+ echo "Subtract $2 ns from SCT RDOs"
+ timestr="--subtractTime $2"
+ shift;
+ shift;;
+
-g | --geom)
geomstr="--geom $2";
shift;
@@ -115,6 +121,7 @@ echo `date` - $HOSTNAME
echo "File: $file_name"
echo "Geom: $geomstr"
echo "Gain: $gainstr"
+echo "Time: $timestr"
echo "Release: $release_directory"
echo "Output: $output_directory"
echo "Starting: $starting_directory"
@@ -178,7 +185,10 @@ cd "$file_stem"
#
# Run job
#
-faser_digi.py $geomstr $gainstr $tagstr "$file_path"
+faser_digi.py $geomstr $gainstr $timestr $tagstr "$file_path"
+digi_code=$?
+echo "Return code: $digi_code"
+
#
# Print out ending time
date
@@ -189,25 +199,52 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-RDO.root to $outdest"
- mkdir -p $outdest
- eos cp *-RDO.root ${outdest}/ || true
+ thefile=`ls *-RDO.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-RDO.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
- mkdir -p $logdest
+ eos mkdir -p $logdest
eos cp $logfile $logdest/$logfile
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
- mkdir -p $outdest
+ eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $digi_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
index 98c34d0eb1f3f5de958e6ab4f8438d4e9405ee54..50278c33e2435236fff3bc2b6a94c13aa5794817 100755
--- a/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faser_digi_merge.sh
@@ -2,10 +2,11 @@
# Used with a condor file to submit to vanilla universe
#
# Usage:
-# submit_faser_digi_merge.sh [--highGain] dirpath slice nfiles [release_directory] [working_directory]
+# submit_faser_digi_merge.sh dirpath slice nfiles [release_directory] [working_directory]
#
# Options:
-# --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
+# --digiTag <tag> - override digitization tag for calo gain
+# --partial - allow missing files in merge
# --geom - geometry setting
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
@@ -28,6 +29,7 @@ SECONDS=0
#
# Job options strings
gainstr=""
+timestr=""
partialstr=""
geomstr=""
#
@@ -41,6 +43,12 @@ do
shift;
shift;;
+ --subtractTime)
+ echo "Subtract $2 ns from SCT RDOs"
+ timestr="--subtractTime $2"
+ shift;
+ shift;;
+
--partial)
echo "Allowing partial merge"
partialstr="--partial"
@@ -146,6 +154,7 @@ echo `date` - $HOSTNAME
echo "Directory: $dir_path"
echo "Geom: $geomstr"
echo "Gain: $gainstr"
+echo "Time: $timestr"
echo "Slice: $slice"
echo "NFiles: $nfiles"
echo "Release: $release_directory"
@@ -212,7 +221,9 @@ cd "$file_stem"
#
# Run job
#
-faser_digi_merge.py $partialstr $geomstr $gainstr $tagstr --slice $slice --files $nfiles $dir_path
+faser_digi_merge.py $partialstr $geomstr $gainstr $timestr $tagstr --slice $slice --files $nfiles $dir_path
+digi_code=$?
+echo "Return code: $digi_code"
#
# Print out ending time
date
@@ -223,25 +234,52 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-RDO.root to $outdest"
- mkdir -p $outdest
- eos cp *-RDO.root ${outdest}/ || true
+ thefile=`ls *-RDO.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-RDO.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
- mkdir -p $logdest
+ eos mkdir -p $logdest
eos cp $logfile $logdest/$logfile
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
- mkdir -p $outdest
+ eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $digi_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index ff8f1d3792a666bf92dca51b16bf7adf58c7500d..055dd7200201ce69823d798609cb991d85f4e794 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -40,16 +40,22 @@ def faser_pgparser():
help="Specify PDG ID of particle (note plus/minus different) or list (e.g.: --pid -13 13)")
parser.add_argument("--mass", default=105.66, type=float,
help="Specify particle mass (in MeV)")
- parser.add_argument("--radius", default=100., type=float,
+ parser.add_argument("--radius", default=100., type=float, nargs="*",
help="Specify radius (in mm)")
parser.add_argument("--angle", default=0.005, type=float_or_none,
help="Specify angular width (in Rad)")
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",
+ parser.add_argument("--hist_name", default="",
help="Specify energy sampling name for hist sampler file.root:hist")
parser.add_argument("--minE", default=10., type=float,
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index cc1628efdb9245e32b7abbbc9ca60a60b2151fad..169badecf7cf2a392a95b098a858634c5a2523c6 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -101,24 +101,50 @@ 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
print(f"Using pid: {args.pid} => {pidarg}")
+ if isinstance(args.radius, list) and len(args.radius) == 1:
+ args.radius = args.radius[0]
+
+
# 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 +158,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/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index f295da6bdc7b914a47b817a3e3d770dc9731e20a..21e1d51a6af7d92ebac09d051e80c919d1c39687 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -315,5 +315,11 @@ b = time.time()
from AthenaCommon.Logging import log
log.info(f"Finish execution in {b-a} seconds")
-sys.exit(int(sc.isFailure()))
+# Signal errors
+if sc.isSuccess():
+ log.info("Execution succeeded")
+ sys.exit(0)
+else:
+ log.info("Execution failed, return 1")
+ sys.exit(1)
diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
index 20b39c9178b01a3e0252b2cb9351c39359939106..f75b144a1bf2d0a725692dbaa82be34ad788cd5c 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
+++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh
@@ -8,7 +8,10 @@
# --out - specify output location (in EOS) to copy output HITS file
# --log - specify output location (in EOS) for log file
# --geom - specify geometry
+#
+# Monte Carlo options:
# --isMC - needed for MC reco
+# --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)
@@ -25,6 +28,9 @@
# Keep track of time
SECONDS=0
#
+# Job option strings
+gainstr=""
+#
# Parse command-line options
while [ -n "$1" ]
do
@@ -45,9 +51,16 @@ do
shift;;
--isMC)
+ echo "Set isMC true"
ismc=1
shift;;
+ --caloTag)
+ echo "Override calo digi tag with $2"
+ gainstr="--MC_calibTag $2"
+ shift;
+ shift;;
+
--) # End of options
shift; # Eat this
break;; # And stop parsing
@@ -119,6 +132,7 @@ echo `date` - $HOSTNAME
echo "File: $file_name"
echo "Filepath: $file_path"
echo "Geom: $geom"
+echo "Gain: $gainstr"
echo "Release: $release_directory"
echo "Output: $output_directory"
echo "Starting: $starting_directory"
@@ -216,7 +230,9 @@ else
mcstr="--isMC"
fi
#
-faser_reco.py "--nevents=$nevents" $geomstr $tagstr $mcstr "$file_path"
+faser_reco.py "--nevents=$nevents" $geomstr $tagstr $mcstr $gainstr "$file_path"
+reco_code=$?
+echo "Return code: $reco_code"
#
# Print out ending time
date
@@ -230,11 +246,27 @@ if ! [ -z "$outdest" ]
then
echo "Output directory:"
ls -l
- echo "copy *-xAOD.root to $outdest"
- eos mkdir -p $outdest
- # Keep this line from stopping script, so we might get a log file
- # || true ensures script continues even if copy fails
- eos cp *-xAOD.root ${outdest}/ || true
+ thefile=`ls *-xAOD.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ # Keep this line from stopping script, so we might get a log file
+ # || true ensures script continues even if copy fails
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that it worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-xAOD.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Copy log file second
@@ -255,3 +287,11 @@ then
eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error is calypso failed
+if [ $reco_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
diff --git a/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh b/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
index 83a2470f1e93bd9703e0e44c5abca3a9aeebe7c1..6e4d4df2eac92d4836c2b77104093ab2dd81c01c 100755
--- a/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
+++ b/Control/CalypsoExample/Simulation/scripts/submit_faser_simulate.sh
@@ -35,6 +35,11 @@ do
xangle=1
shift;; # This 'eats' the argument
+ --fshift)
+ echo "Applying FLUKA crossing-angle shift"
+ xangle=2
+ shift;;
+
-l | --log)
logdest="$2";
shift;
@@ -183,9 +188,16 @@ cd "${file_stem}"
#fi
if [[ -z "$xangle" ]]; then
faser_simulate.py --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
+elif [[ $xangle == 1 ]]; then
+ faser_simulate.py --yangle -0.000160 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
else
- faser_simulate.py --yangle -0.000150 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ echo "Applying double crossing-angle shift for FLUKA!"
+ faser_simulate.py --yangle -0.000310 --yshift 12.0 --skip "$skip_events" -n "$nevts" "$infile" "$outfile"
+ sim_code=$?
fi
+echo "Return code: $sim_code"
#
# Print out ending time
date
@@ -196,16 +208,34 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-HITS.root to $outdest"
- eos mkdir -p $outdest
- eos cp *-HITS.root ${outdest}/ || true
+ thefile=`ls *-HITS.root`
+ if [ $? -eq 0 ]; then
+ echo "copy $thefile to $outdest"
+ eos mkdir -p $outdest
+ eos cp $thefile ${outdest}/${thefile} || true
+
+ # Check that this worked
+ eos ls ${outdest}/${thefile} > /dev/null
+ if [ $? -eq 0 ]; then
+ echo "file $thefile copied to $outdest"
+ copy_code=0
+ else
+ echo "didnt find $thefile in $outdest !"
+ copy_code=1
+ fi
+ else
+ echo "ls *-HITS.root returned nothing!"
+ copy_code=1
+ fi
fi
#
# Also copy log file
if ! [ -z "$logdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $logdest"
eos mkdir -p $logdest
@@ -213,9 +243,17 @@ then
elif ! [ -z "$outdest" ]
then
cd ..
+ echo "Working directory:"
ls -l
echo "copy $logfile to $outdest"
eos mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+# Make sure to return an error is calypso failed
+if [ $sim_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
+
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/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index 6fcd6829b29ea474d5c7380526d3c46bab9c6871..19b1db116991c082d0771699813df2fe94c0a1d9 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -35,12 +35,17 @@ class RadialPosSampler(Sampler):
def r(self):
"r position sampler"
- fwhm = 2*self.radius
- sig = fwhm/(2 * sqrt(2 * log(2)))
-
- if self.radius < 0:
- return sqrt(random.uniform(0, abs(self.radius**2)))
+ #fwhm = 2*self.radius
+ #sig = fwhm/(2 * sqrt(2 * log(2)))
+
+ if isinstance(self.radius, list) and len(self.radius) == 2:
+ # If a list of length 2, generate uniformally over an annulus from r[0] to r[1]
+ return sqrt(random.uniform(abs(self.radius[0]**2), abs(self.radius[1]**2)))
+ elif self.radius < 0:
+ # Generate uniformally up to |r| if r is < 0
+ return sqrt(random.uniform(0, abs(self.radius**2)))
else:
+ # Else generate as a Gaussian with widht r
x = random.gauss(0, self.radius)
y = random.gauss(0, self.radius)
return sqrt(x**2 + y**2)
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 0f725fe9a5a4bf2dc63b915ccf9768b1ffd60c0f..27b30e25df16dd8baf1c4ce4b82b6a9e6a0f7b6f 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -22,7 +22,9 @@ parser.add_argument("path",
parser.add_argument("--slice", type=int, default=0,
help="Specify ordinal output file to produce")
parser.add_argument("--files", type=int, default=1,
- help="Specify files per slice")
+ help="Specify reco files per slice")
+parser.add_argument("--merge", type=int, default=1,
+ help="Specify merged files per reco file (MC only)")
parser.add_argument("--last", type=int, default=0,
help="Specify last file in slice (normally --files)")
@@ -41,7 +43,15 @@ parser.add_argument("--isMC", action='store_true',
help="Running on digitised MC rather than data")
parser.add_argument("--partial", action='store_true',
help="Allow partial input files")
+
+parser.add_argument("--unblind", action='store_true',
+ help="Don't apply signal blinding (default: do)")
+parser.add_argument("--fluka", action='store_true',
+ help="Add FLUKA weights to ntuple")
+parser.add_argument("--genie", action='store_true',
+ help="Add Genie weights to ntuple")
+
args = parser.parse_args()
from pathlib import Path
@@ -53,21 +63,42 @@ filelist = []
# If this is a directory, need to create file list
if filepath.is_dir():
- # Parsing MC is tricky
- if args.isMC:
- print("Monte Carlo not supported yet!")
- sys.exit(0)
-
- # Use expected data pattern to find files
+ # Use expected data pattern to find files (data only)
runstr = filepath.stem
- start = args.slice * args.files
- if args.last > 0:
- end = start + args.last
+
+ # Make list of segments to search for
+ seglist = []
+ if args.merge > 1:
+ start = args.slice * args.files * args.merge
+
+ # Number of files to combine
+ if args.last > 0:
+ num = args.last
+ else:
+ num = args.files
+
+ # Make segment list
+ for i in range(start, start+num*args.merge, args.merge):
+ seg = f"{i:05d}-{(i+args.merge-1):05d}"
+ seglist.append(seg)
+
else:
- end = start + args.files
+ start = args.slice * args.files
+ if args.last > 0:
+ end = start + args.last
+ else:
+ end = start + args.files
+
+ seglist = [f'{seg:05d}' for seg in range(start, end)]
+
+
+ for seg in seglist:
+
+ if args.isMC:
+ searchstr = f"FaserMC-*-{seg}-*xAOD.root"
+ else:
+ searchstr = f"Faser-Physics-{runstr}-{seg}-*xAOD.root"
- for seg in range(start, end):
- searchstr = f"Faser-Physics-{runstr}-{seg:05d}-*xAOD.root"
flist = list(filepath.glob(searchstr))
if len(flist) == 0:
print(f"Didn't find file {searchstr}!")
@@ -86,21 +117,41 @@ 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]
lastfile = Path(filelist[-1])
laststem = str(lastfile.stem)
lastseg = laststem.split('-')[3]
+ 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}")
+ print(f"Blind = {not args.unblind}")
# Find any tags
- tagstr = firststem.replace(f"Faser-Physics-{runstr}-{firstseg}", "")
+ 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"Faser-Physics-{runstr}-{firstseg}-{lastseg}"
+ outfile = f"{firstfaser}-{firstshort}-{runstr}-{firstseg}-{lastseg}"
# This will include the leading -
if len(tagstr) > 0:
@@ -111,6 +162,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:
@@ -161,12 +213,11 @@ Configurable.configurableRun3Behavior = True
# Configure
ConfigFlags.Input.Files = filelist
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ConfigFlags.Input.isMC = args.isMC
if args.isMC:
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions
- ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
else:
ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions
- ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
@@ -184,14 +235,21 @@ acc.merge(PoolReadCfg(ConfigFlags))
# algorithm
from NtupleDumper.NtupleDumperConfig import NtupleDumperAlgCfg
if args.isMC:
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, UseFlukaWeights=True))
+ if args.genie:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, UseGenieWeights=True))
+ elif args.fluka:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, UseFlukaWeights=True))
+ else:
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile))
+
else:
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile))
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind)))
-from AthenaConfiguration.ComponentFactory import CompFactory
-AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
-AthenaEventLoopMgr.EventPrintoutInterval=1000
-acc.addService(AthenaEventLoopMgr)
+if not args.verbose:
+ from AthenaConfiguration.ComponentFactory import CompFactory
+ AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+ AthenaEventLoopMgr.EventPrintoutInterval=1000
+ acc.addService(AthenaEventLoopMgr)
# Hack to avoid problem with our use of MC databases when isMC = False
if not args.isMC:
@@ -217,4 +275,10 @@ b = time.time()
log.info(f"Finish execution in {b-a} seconds")
# Success should be 0
-sys.exit(int(sc.isFailure()))
+if sc.isSuccess():
+ log.info("Execution succeeded")
+ sys.exit(0)
+else:
+ log.info("Execution failed, return 1")
+ sys.exit(1)
+
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh b/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh
index d32998724cb479aeff240dee489257b9f1a0f3fd..9e3273dbb970f5b5d80958f67656490c045a7612 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh
+++ b/PhysicsAnalysis/NtupleDumper/scripts/submit_faser_ntuple_maker.sh
@@ -10,6 +10,9 @@
# --log - specify output location (in EOS) for log file
# --isMC - needed for MC reco
# --partial - allow missing files
+# --merge - merge factor of reco files (for MC)
+# --fluka - create fluka weights
+# --genie - create genie weights
#
# dirpath - full directory path to HITS files
# slice - ordinal output file number
@@ -28,6 +31,10 @@ SECONDS=0
# Defaults
ismc=""
partialstr=""
+mergestr=""
+flukastr=""
+geniestr=""
+unblindstr=""
#
# Parse command-line options
while [ -n "$1" ]
@@ -52,6 +59,23 @@ do
partialstr="--partial"
shift;;
+ --merge)
+ mergestr="--merge $2";
+ shift;
+ shift;;
+
+ --fluka)
+ flukastr="--fluka";
+ shift;;
+
+ --genie)
+ geniestr="--genie";
+ shift;;
+
+ --unblind)
+ unblindstr="--unblind";
+ shift;;
+
--) # End of options
shift; # Eat this
break;; # And stop parsing
@@ -140,6 +164,7 @@ echo `date` - $HOSTNAME
echo "Directory: $dir_path"
echo "Slice: $slice"
echo "NFiles: $nfiles"
+echo "Merge: $mergestr"
echo "Release: $release_directory"
echo "Output: $output_directory"
echo "Starting: $starting_directory"
@@ -219,7 +244,9 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
# Run job
#
-faser_ntuple_maker.py $last_file_str $partialstr $tagstr $ismc --slice $slice --files $nfiles $dir_path
+faser_ntuple_maker.py $last_file_str $partialstr $tagstr $ismc --slice $slice --files $nfiles $mergestr $flukastr $geniestr $unblindstr $dir_path
+ntup_code=$?
+echo "Return code: $ntup_code"
#
# Print out ending time
date
@@ -250,3 +277,6 @@ then
mkdir -p $outdest
eos cp $logfile $outdest/$logfile
fi
+
+# Make sure to return an error if ntuple failed
+exit $ntup_code
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index b61be213ac21f4c7cb568b9bb97177260c266090..3f476edb40480e5871a2df39a3cd2695ea894e4c 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -106,6 +106,7 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_preshowerCalibratedContainer.initialize());
ATH_CHECK(m_ecalCalibratedContainer.initialize());
+ ATH_CHECK(m_eventInfoKey.initialize());
ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
@@ -160,8 +161,8 @@ StatusCode NtupleDumperAlg::initialize()
//WAVEFORMS
addWaveBranches("VetoNu",2,4);
- addWaveBranches("VetoSt1",2,6);
- addWaveBranches("VetoSt2",1,14);
+ addWaveBranches("VetoSt1",1,14);
+ addWaveBranches("VetoSt2",2,6);
addWaveBranches("Timing",4,8);
addWaveBranches("Preshower",2,12);
addWaveBranches("Calo",4,0);
@@ -197,6 +198,7 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("TrackSegment_pz", &m_trackseg_pz);
//TrackCollection
+ m_tree->Branch("Track_PropagationError", &m_propagationError, "Track_PropagationError/I");
m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
m_tree->Branch("Track_Chi2", &m_Chi2);
m_tree->Branch("Track_nDoF", &m_DoF);
@@ -304,6 +306,24 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I");
m_tree->Branch("CrossSection", &m_crossSection, "crossSection/D");
+
+ // first 10 truth particles
+
+ m_tree->Branch("truth_P", &m_truth_P);
+ m_tree->Branch("truth_px", &m_truth_px);
+ m_tree->Branch("truth_py", &m_truth_py);
+ m_tree->Branch("truth_pz", &m_truth_pz);
+ m_tree->Branch("truth_m", &m_truth_m);
+ m_tree->Branch("truth_pdg", &m_truth_pdg);
+
+ m_tree->Branch("truth_prod_x", &m_truth_prod_x);
+ m_tree->Branch("truth_prod_y", &m_truth_prod_y);
+ m_tree->Branch("truth_prod_z", &m_truth_prod_z);
+
+ m_tree->Branch("truth_dec_x", &m_truth_dec_x);
+ m_tree->Branch("truth_dec_y", &m_truth_dec_y);
+ m_tree->Branch("truth_dec_z", &m_truth_dec_z);
+
// for mother + daughter particle truth infomation
m_tree->Branch("truthM_P", &m_truthM_P);
@@ -410,6 +430,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
bool trig_coincidence_vetoesORpreshower_and_calo = false;
if ( ((m_tap&2)!=0) && ((m_tap&1)!=0) ) trig_coincidence_vetoesORpreshower_and_calo = true;
+ bool trig_calo = (m_tap & 1);
+
// for random trigger, store charge of scintillators in histograms
if (trig_random) {
// Read in Waveform containers
@@ -462,7 +484,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
return StatusCode::SUCCESS; // finished with this randomly triiggered event
- } else if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_coincidence_timing_and_calo || trig_coincidence_vetoesORpreshower_and_calo) ) {
+ // Try adding calo-only trigger
+ } else if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_coincidence_timing_and_calo || trig_coincidence_vetoesORpreshower_and_calo || trig_calo) ) {
// don't process events that fail to activate coincidence triggers
return StatusCode::SUCCESS;
}
@@ -546,30 +569,64 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
{
- for (auto particle : *truthParticleContainer)
- {
- if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 )
- {
-
- if ( particle->pdgId() == 32) // mother particle (A')
- {
- m_truthM_P.push_back(particle->p4().P());
- m_truthM_px.push_back(particle->p4().X());
- m_truthM_py.push_back(particle->p4().Y());
- m_truthM_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthM_x.push_back(particle->prodVtx()->x());
- m_truthM_y.push_back(particle->prodVtx()->y());
- m_truthM_z.push_back(particle->prodVtx()->z());
- } else {
+ int ipart(0);
+ for (auto particle : *truthParticleContainer)
+ {
+
+ // loop over first 10 truth particles (for non A' samples)
+
+ if (ipart++ < 10) {
+
+ m_truth_P.push_back(particle->p4().P());
+ m_truth_px.push_back(particle->p4().X());
+ m_truth_py.push_back(particle->p4().Y());
+ m_truth_pz.push_back(particle->p4().Z());
+ m_truth_m.push_back(particle->m());
+ m_truth_pdg.push_back(particle->pdgId());
+
+ if ( particle->hasProdVtx()) {
+ m_truth_prod_x.push_back(particle->prodVtx()->x());
+ m_truth_prod_y.push_back(particle->prodVtx()->y());
+ m_truth_prod_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truth_prod_x.push_back(NaN);
+ m_truth_prod_y.push_back(NaN);
+ m_truth_prod_z.push_back(NaN);
+ }
+
+ if ( particle->hasDecayVtx()) {
+ m_truth_dec_x.push_back(particle->decayVtx()->x());
+ m_truth_dec_y.push_back(particle->decayVtx()->y());
+ m_truth_dec_z.push_back(particle->decayVtx()->z());
+ } else {
+ m_truth_dec_x.push_back(NaN);
+ m_truth_dec_y.push_back(NaN);
+ m_truth_dec_z.push_back(NaN);
+ }
+ }
+
+
+ if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
+
+ if ( particle->pdgId() == 32) { // mother particle (A')
+
+ m_truthM_P.push_back(particle->p4().P());
+ m_truthM_px.push_back(particle->p4().X());
+ m_truthM_py.push_back(particle->p4().Y());
+ m_truthM_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasDecayVtx()) { // decay vertex for A' particle
+ m_truthM_x.push_back(particle->decayVtx()->x());
+ m_truthM_y.push_back(particle->decayVtx()->y());
+ m_truthM_z.push_back(particle->decayVtx()->z());
+ } else {
m_truthM_x.push_back(NaN);
m_truthM_y.push_back(NaN);
m_truthM_z.push_back(NaN);
}
+ }
- }
if ( particle->pdgId() == 11) // daughter particle (positron)
{
m_truthd0_P.push_back(particle->p4().P());
@@ -604,12 +661,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_truthd1_z.push_back(NaN);
}
}
- }
+ }
}
- }
+ }
}
// load in calibrated calo container
+ if (m_doCalib) {
SG::ReadHandle<xAOD::CalorimeterHitContainer> ecalCalibratedContainer { m_ecalCalibratedContainer, ctx };
ATH_CHECK(ecalCalibratedContainer.isValid());
for (auto hit : *ecalCalibratedContainer) {
@@ -649,6 +707,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 };
@@ -673,14 +732,23 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
FillWaveBranches(*ecalContainer);
// enforce blinding such that events that miss the veto layers and have a large calo signal are skipped and not in the output root file
- if ( (!isMC) && m_doBlinding) {
- if ( m_calo_total_E_EM/1000.0 > 10.0 ) { // energy is in MeV so divide by 1000 to compare to 10 GeV
+ // Only blind colliding BCIDs (+/- 1)
+ if ( (!isMC) && m_doBlinding && abs(m_distanceToCollidingBCID) <= 1) {
+ if (m_calo_total_E_EM > m_blindingCaloThreshold ) { // energy is in MeV
if (m_wave_status[4] == 1 and m_wave_status[5] == 1 and m_wave_status[6] == 1 and m_wave_status[7] == 1 and m_wave_status[14] == 1) { // hit status == 1 means it is below threshold. channles 4 and 5 are vetoNu, channels 6,7, and 14 are veto
return StatusCode::SUCCESS;
}
}
}
+ SG::ReadDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx);
+ if (eventInfo->errorState(xAOD::EventInfo_v1::SCT) == xAOD::EventInfo::Error) {
+ m_propagationError = 1;
+ ATH_MSG_DEBUG("NtupleDumper: xAOD::EventInfo::SCT::Error");
+ } else {
+ m_propagationError = 0;
+ }
+
// get geometry context
FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
auto gctx = faserGeometryContext.context();
@@ -766,8 +834,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;
@@ -1219,7 +1295,20 @@ NtupleDumperAlg::clearTree() const
m_truthBarcode = 0;
m_truthPdg = 0;
- m_truthM_P.clear();
+ m_truth_P.clear();
+ m_truth_px.clear();
+ m_truth_py.clear();
+ m_truth_pz.clear();
+ m_truth_m.clear();
+ m_truth_pdg.clear();
+ m_truth_prod_x.clear();
+ m_truth_prod_y.clear();
+ m_truth_prod_z.clear();
+ m_truth_dec_x.clear();
+ m_truth_dec_y.clear();
+ m_truth_dec_z.clear();
+
+ m_truthM_P.clear();
m_truthM_px.clear();
m_truthM_py.clear();
m_truthM_pz.clear();
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 5cc342021993ca07354e6a4924f7e5cfb07ae552..6d4ccc381a4b65831678a657c035ecc4e79c2b65 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -21,6 +21,8 @@
#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "StoreGate/ReadDecorHandle.h"
#include <vector>
@@ -84,7 +86,8 @@ private:
SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
SG::ReadHandleKey<xAOD::WaveformClock> m_ClockWaveformContainer { this, "WaveformClockKey", "WaveformClock", "ReadHandleKey for ClockWaveforms Container"};
- ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+ SG::ReadDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
@@ -98,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" };
@@ -106,6 +111,8 @@ private:
DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
+ DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"};
+
double m_baseEventCrossSection {1.0};
const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -189,6 +196,7 @@ private:
mutable std::vector<double> m_trackseg_pz;
mutable int m_longTracks;
+ mutable int m_propagationError;
mutable std::vector<double> m_Chi2;
mutable std::vector<double> m_DoF;
mutable std::vector<double> m_xup;
@@ -274,7 +282,7 @@ private:
mutable std::vector<double> m_truthM_py;
mutable std::vector<double> m_truthM_pz;
- mutable std::vector<double> m_truthM_x;
+ mutable std::vector<double> m_truthM_x; // decay vertex of A'
mutable std::vector<double> m_truthM_y;
mutable std::vector<double> m_truthM_z;
@@ -284,7 +292,7 @@ private:
mutable std::vector<double> m_truthd0_py;
mutable std::vector<double> m_truthd0_pz;
- mutable std::vector<double> m_truthd0_x;
+ mutable std::vector<double> m_truthd0_x; // production vertex for daughter particles
mutable std::vector<double> m_truthd0_y;
mutable std::vector<double> m_truthd0_z;
@@ -298,6 +306,27 @@ private:
mutable std::vector<double> m_truthd1_y;
mutable std::vector<double> m_truthd1_z;
+ // first 10 truth particles
+
+ mutable std::vector<double> m_truth_P;
+ mutable std::vector<double> m_truth_px;
+ mutable std::vector<double> m_truth_py;
+ mutable std::vector<double> m_truth_pz;
+ mutable std::vector<double> m_truth_m;
+
+ mutable std::vector<double> m_truth_dec_x; // components of decay vertex (mm)
+ mutable std::vector<double> m_truth_dec_y;
+ mutable std::vector<double> m_truth_dec_z;
+
+ mutable std::vector<double> m_truth_prod_x; // components of production vertex (mm)
+ mutable std::vector<double> m_truth_prod_y;
+ mutable std::vector<double> m_truth_prod_z;
+
+ mutable std::vector<int> m_truth_pdg; // pdg of first 10 truth particles
+
+
+
+
mutable double m_truthLeptonMomentum;
mutable int m_truthBarcode;
mutable int m_truthPdg;
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index ea03a77f16d12cd4cb45c99c96a20c1fbc8323e7..1869860e84a401fffbfe6a588d6b8e78ba702e9a 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -3,6 +3,7 @@
#include "ScintSimEvent/ScintHitIdHelper.h"
#include <map>
+#include <cmath>
ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name,
@@ -110,42 +111,68 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
waveforms = m_digiTool->create_waveform_map(m_preshowerID);
}
+ // Sum energy for each channel (i.e. identifier)
+ std::map<Identifier, float> esum;
+
+ Identifier id;
+ for (const auto& hit : *scintHitHandle) {
+ if (first.isTrigger()) {
+ Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+ // Need to do something better here, but for now just fill both
+ id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+ esum[id] += hit.energyLoss();
+ id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+ esum[id] += hit.energyLoss();
+ } else {
+ // All others have only 1 PMT
+ // Use detector id (not hit id) throughout
+ id = hit.getIdentifier();
+ esum[id] += hit.energyLoss();
+ }
+ }
+
// Loop over time samples
- // Should really do sums first, as repeating these in the loop is a waste
for (const auto& tk : m_timekernel) {
std::map<Identifier, float> counts;
- // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
- Identifier id;
- for (const auto& hit : *scintHitHandle) {
-
- // Special handling for trigger scintillator
- if (first.isTrigger()) {
- // Hits are created per plate
- // Need to split between 2 PMTs
-
- // Identifier from hit is plate ID
- Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
-
- // Need to do something better here, but for now just fill both
- id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
- counts[id] += tk * hit.energyLoss();
- id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
- counts[id] += tk * hit.energyLoss();
-
- } else {
- // All others have only 1 PMT
- // Use detector id (not hit id) throughout
- id = hit.getIdentifier();
- counts[id] += tk * hit.energyLoss();
- }
+// // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
+// Identifier id;
+// for (const auto& hit : *scintHitHandle) {
+//
+// // Special handling for trigger scintillator
+// if (first.isTrigger()) {
+// // Hits are created per plate
+// // Need to split between 2 PMTs
+//
+// // Identifier from hit is plate ID
+// Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+//
+// // Need to do something better here, but for now just fill both
+// id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+// counts[id] += tk * hit.energyLoss();
+// id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+// counts[id] += tk * hit.energyLoss();
+//
+// } else {
+// // All others have only 1 PMT
+// // Use detector id (not hit id) throughout
+// id = hit.getIdentifier();
+// counts[id] += tk * hit.energyLoss();
+// }
+// }
+
+ // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+ for (const auto& e : esum) {
+ // Convert hit id to Identifier and store
+ id = e.first;
+ counts[id] = tk * e.second;
}
// Subtract count from basleine and add result to correct waveform vector
for (const auto& c : counts) {
- unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = baseline - c.second;
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = std::round(baseline - c.second);
if (value < 0) {
ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.cxx b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.cxx
index ea353da6c61a80535c20a3f92f3598a227b65b0f..bada70e4db0f3e105788e0042f3986ffef3e4bef 100644
--- a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.cxx
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.cxx
@@ -170,6 +170,20 @@ StatusCode FaserSCT_SurfaceChargesGenerator::initialize() {
// ----------------------------------------------------------------------
StatusCode FaserSCT_SurfaceChargesGenerator::finalize() {
ATH_MSG_DEBUG("FaserSCT_SurfaceChargesGenerator::finalize()");
+
+ if (m_tofNum > 0) {
+ float mean_TOF = m_tofSum / m_tofNum;
+ float rms_TOF = sqrt(m_tofSum2 / m_tofNum - mean_TOF*mean_TOF);
+ ATH_MSG_INFO("TOF hits: " << m_tofNum);
+ ATH_MSG_INFO("Mean TOF: " << mean_TOF);
+ ATH_MSG_INFO("RMS TOF: " << rms_TOF);
+
+ }
+
+ if (m_tfix > -998.) {
+ ATH_MSG_INFO("Used fixed value: " << m_tfix.value());
+ }
+
return StatusCode::SUCCESS;
}
@@ -335,13 +349,22 @@ void FaserSCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* ele
float timeOfFlight{p_eventTime + hitTime(phit)};
// Kondo 19/09/2007: Use the coordinate of the center of the module to calculate the time of flight
- timeOfFlight -= (element->center().mag()) / CLHEP::c_light;
+ // Torrence 19/02/2023: Use global z coordinate rather than distance to origin for FASER
+ timeOfFlight -= (element->center().z()) / CLHEP::c_light;
// !< extract the distance to the origin of the module to Time of flight
// !< timing set from jo to adjust (subtract) the timing
if (m_tsubtract > -998.) {
timeOfFlight -= m_tsubtract;
}
+
+ ATH_MSG_DEBUG("Time of flight: " << timeOfFlight);
+
+ // Keep some stats
+ m_tofNum += 1;
+ m_tofSum += timeOfFlight;
+ m_tofSum2 += (timeOfFlight * timeOfFlight);
+
// ---**************************************
const CLHEP::Hep3Vector pos{phit.localStartPosition()};
@@ -512,6 +535,7 @@ void FaserSCT_SurfaceChargesGenerator::processSiHit(const SiDetectorElement* ele
const float sdist{static_cast<float>(design->scaledDistanceToNearestDiode(position))}; // !< dist on the surface from the hit point to the nearest strip (diode)
const float t_surf{surfaceDriftTime(2.0 * sdist)}; // !< Surface drift time
const float totaltime{(m_tfix > -998.) ? m_tfix.value() : t_drift + timeOfFlight + t_surf}; // !< Total drift time
+ ATH_MSG_VERBOSE(std::fixed << std::setprecision(4) << "Surface time: " << t_surf << " Drift time: " << t_drift << " TOF: " << timeOfFlight << " Total: " << totaltime);
inserter(SiSurfaceCharge(position, SiCharge(q1, totaltime, hitproc, trklink)));
} else {
ATH_MSG_VERBOSE(std::fixed << std::setprecision(8) << "Local position (phi, eta, depth): ("
diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.h b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.h
index b84f584c52ad1ac182f237d62a7639aeb0fe2f13..4e9ec5775edc2c3b90fd095c411e66435bc7d754 100644
--- a/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.h
+++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/src/FaserSCT_SurfaceChargesGenerator.h
@@ -137,6 +137,12 @@ class FaserSCT_SurfaceChargesGenerator : public extends<AthAlgTool, ISCT_Surface
bool m_SurfaceDriftFlag{false}; //!< surface drift ON/OFF
+ // Keep track of TOF
+ // These are updated in processSiHit which is const, so make mutable
+ mutable int m_tofNum{0};
+ mutable float m_tofSum{0.};
+ mutable float m_tofSum2{0.};
+
// -- Histograms
// TProfile* m_h_efieldz{nullptr};
// TH1F* m_h_efield{nullptr};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
index 0c890e1eb0542037393884a0d181abe9799c6664..5d6ed5c4f3170cfc57b058d6bc9432329c1da561 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
+++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py
@@ -107,9 +107,7 @@ def CKF2Cfg(flags, **kwargs):
ckf.PerformanceWriter = False
ckf.nMax = 10
- # Use larger chi2 cut until alignment improves
- # ckf.chi2Max = 25
- ckf.chi2Max = 100000
+ ckf.chi2Max = 25
# Protect against running out of memory on busy events
ckf.maxSteps = 5000
acc.addEventAlgo(ckf)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 3647ac98abdd9e5176fb0902fb368e89dc741dcb..f990be9458e0ee2ac0a5bdef8b670703f1424c26 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -2,6 +2,7 @@
#include "StoreGate/ReadHandle.h"
#include "StoreGate/ReadCondHandleKey.h"
+#include "StoreGate/WriteDecorHandle.h"
#include "TrackerSpacePoint/FaserSCT_SpacePointCollection.h"
#include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
#include "TrackerIdentifier/FaserSCT_ID.h"
@@ -28,9 +29,9 @@
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
-
-
#include "Acts/EventData/Measurement.hpp"
+#include "Acts/Propagator/PropagatorError.hpp"
+#include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp"
size_t CKF2::TrajectoryInfo::nClusters {0};
@@ -49,6 +50,7 @@ StatusCode CKF2::initialize() {
ATH_CHECK(m_kalmanFitterTool1.retrieve());
ATH_CHECK(m_createTrkTrackTool.retrieve());
ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_eventInfoKey.initialize());
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool.retrieve());
}
@@ -78,6 +80,8 @@ StatusCode CKF2::execute() {
const EventContext& ctx = Gaudi::Hive::currentContext();
m_numberOfEvents++;
+ SG::WriteDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx);
+
SG::WriteHandle trackContainer{m_trackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
@@ -150,6 +154,14 @@ StatusCode CKF2::execute() {
std::list<TrajectoryInfo> allTrajectories;
for (auto &result : results) {
if (not result.ok()) {
+ // TODO use status bits for different errors
+ // result.error() == Acts::CombinatorialKalmanFilterError::NoTrackFound
+ if (result.error() == Acts::PropagatorError::StepCountLimitReached ||
+ result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) {
+ if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) {
+ ATH_MSG_WARNING (" cannot set error state for SCT");
+ }
+ }
continue;
}
CKFResult ckfResult = result.value();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 2111afcf6b52d46b45e1e2539d105290b937689e..170f0353904bef82095146c50bf98d5721ef2138 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -21,6 +21,7 @@
#include "KalmanFitterTool.h"
#include <boost/dynamic_bitset.hpp>
#include "CreateTrkTrackTool.h"
+#include "xAODEventInfo/EventInfo.h"
using ConstTrackStateProxy = Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>;
using ClusterSet = boost::dynamic_bitset<>;
@@ -148,6 +149,7 @@ private:
ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
Gaudi::Property<bool> m_isMC {this, "isMC", false};
SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "CKFTrackCollection", "Output track collection name" };
+ SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
};
#endif // FASERACTSKALMANFILTER_CKF2_H
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index c30351902a1b79e9a58ad33d5f58a5d04ed13faa..37cb7af183c2d89b9ada3579f76820bfe6487c5c 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -51,7 +51,7 @@ public:
virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0;
/// Generate random baseline
- virtual unsigned int generate_baseline(int mean, int rms) const = 0;
+ virtual float generate_baseline(float mean, float rms) const = 0;
/// Create structure to store pulse for each channel
template <class T>
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index c62d7f8a753490907bd1e8a513cad5341a2a46be..3a48027d57e3710a2ec592aed6a9ed749b9f6848 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -22,7 +22,7 @@ WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
m_nsamples = 600;
- m_random = new TRandom3();
+ m_random = new TRandom3(0);
return StatusCode::SUCCESS;
}
@@ -41,7 +41,7 @@ WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
return timekernel;
}
-unsigned int
-WaveformDigitisationTool::generate_baseline(int mean, int rms) const {
+float
+WaveformDigitisationTool::generate_baseline(float mean, float rms) const {
return m_random->Gaus(mean, rms);
}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index e2dd5169152845824927baeeae7ce8fc36ab46f8..b6554c337830fbd2183752cc1b05123ea827b4f1 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -32,7 +32,7 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
std::vector<float> evaluate_timekernel(TF1* kernel) const;
/// Generate random baseline
- unsigned int generate_baseline(int mean, int rms) const;
+ float generate_baseline(float mean, float rms) const;
private: