diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 055dd7200201ce69823d798609cb991d85f4e794..03817efb8d9f8852a1dda650a78f7f4de8283d4e 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -44,6 +44,10 @@ def faser_pgparser():
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("--xpos", default=None, type=float,
+ help="Specify x position of particles (in mm)")
+ parser.add_argument("--ypos", default=None, type=float,
+ help="Specify y position of particles (in mm)")
parser.add_argument("--zpos", default=None, type=float,
help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index 169badecf7cf2a392a95b098a858634c5a2523c6..7ec05e2aaf4df3d434577b3d7f423bcf56856722 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -148,9 +148,15 @@ if __name__ == '__main__':
# -1000 is safely upstream of detector (to be checked)
# Note zpos is in mm!
- if args.zpos:
+ if args.zpos != None:
sg_dict["z"] = args.zpos
+ if args.xpos != None:
+ sg_dict["x"] = args.xpos
+
+ if args.ypos != None:
+ sg_dict["y"] = args.ypos
+
# Determine energy sampling
if args.sampler == "lin":
sg_dict["energy"] = PG.UniformSampler(args.minE*GeV, args.maxE*GeV)
@@ -159,11 +165,28 @@ if __name__ == '__main__':
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)
+ nargs = len(args.hist_name.split(":"))
+ if nargs == 2:
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+ elif nargs == 3:
+ fname, hname, scale = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH1Sampler(fname, hname, scale)
+ else:
+ print(f"Can't parse histogram {args.hist_name}!")
+ sys.exit(1)
+
elif args.sampler == "hist2D":
- fname, hname = args.hist_name.split(":")
- sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+ nargs = len(args.hist_name.split(":"))
+ if nargs == 2:
+ fname, hname = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+ elif nargs == 4:
+ fname, hname, scalex, scaley = args.hist_name.split(":")
+ sg_dict["energy"] = PG.TH2Sampler(fname, hname, scalex, scaley)
+ else:
+ print(f"Can't parse histogram {args.hist_name}!")
+ sys.exit(1)
else:
print(f"Sampler {args.sampler} not known!")
sys.exit(1)
@@ -255,8 +278,12 @@ if __name__ == '__main__':
b = time.time()
log.info("Run G4FaserAlg in " + str(b-a) + " seconds")
-#
-# Success should be 0
-#
- sys.exit(not sc.isSuccess())
+
+# 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/Generation/scripts/submit_faser_particlegun.sh b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
index 9b8676752be083d66acaaf0d0fa4744be43050b5..441ef2923b397274b74f53b8530adf592dc3e76f 100755
--- a/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
+++ b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
@@ -185,9 +185,12 @@ cd "${config_file_stem}-${seg_str}"
# Run job
if [[ -z "$tag" ]]; then
faser_particlegun.py "--conf=$config_path" "--segment=$seg_str"
+ gen_code=$?
else
faser_particlegun.py "--conf=$config_path" "--segment=$seg_str" "--tag=$tag"
+ gen_code=$?
fi
+echo "Return code: $gen_code"
#
# Print out ending time
date
@@ -198,25 +201,51 @@ export EOS_MGM_URL=root://eospublic.cern.ch
#
if ! [ -z "$outdest" ]
then
+ echo "Output directory:"
ls -l
- echo "copy *-HITS.root to $outdest"
- 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 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
#
# 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 ..
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 [ $gen_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+ exit 1
+else
+ exit 0
+fi
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 21e1d51a6af7d92ebac09d051e80c919d1c39687..744aafa0b05213073086939fbf2d75138479b6a6 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -29,6 +29,8 @@ parser.add_argument("-r", "--reco", default="",
help="Specify reco tag (to append to output filename)")
parser.add_argument("-n", "--nevents", type=int, default=-1,
help="Specify number of events to process (default: all)")
+parser.add_argument("--skip", type=int, default=0,
+ help="Specify number of events to skip (default: none)")
parser.add_argument("-v", "--verbose", action='store_true',
help="Turn on DEBUG output")
parser.add_argument("--isMC", action='store_true',
@@ -37,6 +39,8 @@ parser.add_argument("--MC_calibTag", default="",
help="Specify tag used to reconstruct MC calo energy: (WAVE-Calibration-01-LG-nofilt, WAVE-Calibration-01-LG, WAVE-Calibration-01-HG-nofilt, or WAVE-Calibration-01-HG) ")
parser.add_argument("--testBeam", action='store_true',
help="Set geometry for 2021 test beam")
+parser.add_argument("--isOverlay", action='store_true',
+ help="Set overlaid data input")
args = parser.parse_args()
@@ -81,6 +85,8 @@ else:
print(f"Starting reconstruction of {filepath.name} with type {runtype}")
if args.nevents > 0:
print(f"Reconstructing {args.nevents} events by command-line option")
+if args.skip > 0:
+ print(f"Skipping {args.skip} events by command-line option")
# Start reconstruction
@@ -103,6 +109,8 @@ else:
ConfigFlags.Input.ProjectName = "data20"
ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Exec.SkipEvents = args.skip
+
# Flags for later
useCKF = True
useCal = False
@@ -173,7 +181,7 @@ acc.merge(PoolWriteCfg(ConfigFlags))
#
# Set up RAW data access
-if args.isMC:
+if args.isMC or args.isOverlay:
from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
acc.merge(PoolReadCfg(ConfigFlags))
else:
@@ -185,22 +193,24 @@ else:
from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
acc.merge(FaserGeometryCfg(ConfigFlags))
-if useLHC:
+if useLHC and not args.isOverlay:
from LHCDataAlgs.LHCDataAlgConfig import LHCDataAlgCfg
acc.merge(LHCDataAlgCfg(ConfigFlags))
# Set up algorithms
-from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
-acc.merge(WaveformReconstructionCfg(ConfigFlags))
+if not args.isOverlay:
+ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+ acc.merge(WaveformReconstructionCfg(ConfigFlags))
-# Calorimeter Energy reconstruction
-if useCal:
- from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
- acc.merge(CalorimeterReconstructionCfg(ConfigFlags, MC_calibTag=args.MC_calibTag))
+ # Calorimeter Energy reconstruction
+ if useCal:
+ from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
+ acc.merge(CalorimeterReconstructionCfg(ConfigFlags, MC_calibTag=args.MC_calibTag))
# Tracker clusters
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
-acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_RDOs"))
+# acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="Pos_SCT_RDOs"))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_RDOs", checkBadChannels=True))
#
# SpacePoints
@@ -226,7 +236,8 @@ if useCKF:
#
# Kalman Filter for tracking
from FaserActsKalmanFilter.CKF2Config import CKF2Cfg
- acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
+ if not args.isOverlay:
+ acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
# Add tracking collection with no IFT
acc.merge(CKF2Cfg(ConfigFlags, maskedLayers=[0, 1, 2], name="CKF_woIFT",
@@ -246,10 +257,10 @@ itemList = [ "xAOD::EventInfo#*"
, "TrackCollection#*"
]
#
-if useLHC:
+if useLHC and not args.isOverlay:
itemList.extend( ["xAOD::FaserLHCData#*", "xAOD::FaserLHCDataAux#*"] )
-if args.isMC:
+if args.isMC and not args.isOverlay:
# Make xAOD versions of truth
from Reconstruction.xAODTruthCnvAlgConfig import xAODTruthCnvAlgCfg
acc.merge(xAODTruthCnvAlgCfg(ConfigFlags))
@@ -266,13 +277,14 @@ tagBuilder = CompFactory.EventInfoTagBuilder()
tagBuilder.PropagateInput=False
acc.addEventAlgo(tagBuilder)
-# Waveform reconstruction output
-from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
-acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
+if not args.isOverlay:
+ # Waveform reconstruction output
+ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
+ acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
-# Calorimeter reconstruction output
-from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
-acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
+ # Calorimeter reconstruction output
+ from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
+ acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
# Check what we have
print( "Writing out xAOD objects:" )
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
index 9f541da7c972ac01170b8acfafc8001b6eea5a14..10549ef0b250b999c7ac66063bdb1d5ae7c90147 100644
--- a/Generators/ParticleGun/python/histsampling.py
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -14,6 +14,7 @@ def load_hist(*args):
"""
Load a histogram from a filename/TFile and histo name. If a single arg is
provided, it has to be a histo object and will be cloned before return.
+
"""
h = None
if len(args) == 1 and issubclass(type(args[0]), ROOT.TH1):
@@ -25,6 +26,7 @@ def load_hist(*args):
#f.Close()
elif type(args[0]) is ROOT.TFile and type(args[1]) is str:
h = args[0].Get(args[1]).Clone()
+
if h is None:
raise Exception("Error in histogram loading from " + args)
else:
@@ -77,7 +79,7 @@ def get_random_bin(globalbins, cheights):
raise Exception("Sample fell outside range of cumulative distribution?!?!")
-def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
+def get_random_x(h, globalbins, cheights, globalbin_to_axisbin, scale=1.):
"""
Choose a random bin via get_random_bin, then pick a uniform random x
point in that bin (without any attempt at estimating the in-bin distribution).
@@ -86,10 +88,10 @@ def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
axisids = globalbin_to_axisbin.get(irand)
assert axisids is not None
xrand = random.uniform(h.GetXaxis().GetBinLowEdge(axisids[0]), h.GetXaxis().GetBinUpEdge(axisids[0]))
- return xrand
+ return scale * xrand
-def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
+def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin, xscale=1., yscale=1.):
"""
Choose a random bin via get_random_bin, then pick a uniform random x,y
point in that bin (without any attempt at estimating the in-bin distribution).
@@ -99,34 +101,71 @@ def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
assert axisids is not None
xrand = random.uniform(h2.GetXaxis().GetBinLowEdge(axisids[0]), h2.GetXaxis().GetBinUpEdge(axisids[0]))
yrand = random.uniform(h2.GetYaxis().GetBinLowEdge(axisids[1]), h2.GetYaxis().GetBinUpEdge(axisids[1]))
- return xrand, yrand
+ return (xscale*xrand), (yscale*yrand)
class TH1(object):
"Minimal wrapper for ROOT TH1, for sampling consistency and easy loading"
def __init__(self, *args):
- self.th1 = load_hist(*args)
+ """ Args can be variable length, but these are now kwargs, so order matters
+
+ fname - filename
+ hname - histogram name
+ xscale - scaling factor for x axis variable
+ """
+ if len(args) > 2:
+ self.th1 = load_hist(*(args[0:2]))
+ else:
+ self.th1 = load_hist(*args)
+
self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+ if len(args) >= 3:
+ self.xscale = float(args[2])
+ else:
+ self.xscale = 1.
+
def GetRandom(self):
"A GetRandom that works for TH1s and uses Python random numbers"
if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th1)
- return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+ return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale)
class TH2(object):
"Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
def __init__(self, *args):
- self.th2 = load_hist(*args)
+ """ Args can be variable length, but these are now kwargs, so order matters
+
+ fname - filename
+ hname - histogram name
+ xscale - scaling factor for x axis variable
+ yscale - scaling factor for y axis variable
+ """
+ if len(args) > 2:
+ self.th2 = load_hist(*(args[0:2]))
+ else:
+ self.th2 = load_hist(*args)
+
self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
+ if len(args) >= 3:
+ self.xscale = float(args[2])
+ else:
+ self.xscale = 1.
+
+ if len(args) >= 4:
+ self.yscale = float(args[3])
+ else:
+ self.yscale = 1.
+
+
def GetRandom(self):
"A GetRandom that works for TH2s"
if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th2)
- return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+ return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale, self.yscale)
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index d6e540ea6cf7a747d9152202cd488c135fa27099..9126b607ca7246480e4129669f3933e236c2c126 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -633,14 +633,23 @@ class EThetaMPhiSampler(MomSampler):
pt = p sin(theta)
"""
- if self._theta is None:
- e,theta = self.energy()
- else:
- e = self.energy()
- theta = self.theta()
-
m = self.mass()
- p = math.sqrt( e**2 - m**2 )
+
+ count = 0
+ e = -1
+ while (e < m and count < 5):
+ count += 1
+ if self._theta is None:
+ e,theta = self.energy()
+ else:
+ e = self.energy()
+ theta = self.theta()
+
+ try:
+ p = math.sqrt( e**2 - m**2 )
+ except Exception:
+ raise Exception(f"Error generating E: {e} m: {m}!")
+
pz = p * math.cos(theta)
pt = p * math.sin(theta)
phi = self.phi()
diff --git a/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx b/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx
index de22698589abfbbbfc2ea9aa6d54feacd4661610..1322903932ee6ffdc8d41f31e3a850d680f1c391 100644
--- a/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx
+++ b/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx
@@ -55,6 +55,10 @@ LHCDataAlg::execute(const EventContext& ctx) const {
lhcDataHandle->set_numBunchBeam2(m_lhcTool->getBeam2Bunches(ctx));
lhcDataHandle->set_numBunchColliding(m_lhcTool->getCollidingBunches(ctx));
+ ATH_MSG_DEBUG("FaserLHCData B1: " << m_lhcTool->getBeam1Bunches(ctx)
+ << " B2: " << m_lhcTool->getBeam2Bunches(ctx)
+ << " Coll: " << m_lhcTool->getCollidingBunches(ctx));
+
// Fill BCID information
// Get the BCID mask
@@ -64,37 +68,68 @@ LHCDataAlg::execute(const EventContext& ctx) const {
SG::ReadHandle<xAOD::EventInfo> xevt(m_eventInfo, ctx);
unsigned int bcid = xevt->bcid();
- int nearest = findNearest(bcid, bcid_mask, 3); // Colliding beams
+ int nearest;
+ if (m_lhcTool->getCollidingBunches(ctx) == 0) {
+ ATH_MSG_INFO("No colliding bunches, can't set nearest");
+ nearest = -3564;
+ } else {
+ nearest = findNearest(bcid, bcid_mask, 3); // Colliding beams
+ }
lhcDataHandle->set_distanceToCollidingBCID(nearest);
ATH_MSG_DEBUG("Found distance of " << nearest << " from BCID " << bcid
<< " to the nearest colliding BCID ");
- nearest = findNearest(bcid, bcid_mask, 1); // Beam1 unpaired
+ if (m_lhcTool->getBeam1Bunches(ctx) == 0) {
+ ATH_MSG_INFO("No beam 1 bunches, can't set nearest");
+ nearest = -3564;
+ } else {
+ nearest = findNearest(bcid, bcid_mask, 1); // Beam1 unpaired
+ }
lhcDataHandle->set_distanceToUnpairedB1(nearest);
ATH_MSG_DEBUG("Found distance of " << nearest << " from BCID " << bcid
<< " to the nearest unpaired B1 ");
- nearest = findNearest(bcid, bcid_mask, 2); // Beam2 unpaired
+ if (m_lhcTool->getBeam2Bunches(ctx) == 0) {
+ ATH_MSG_INFO("No beam 2 bunches, can't set nearest");
+ nearest = -3564;
+ } else {
+ nearest = findNearest(bcid, bcid_mask, 2); // Beam2 unpaired
+ }
lhcDataHandle->set_distanceToUnpairedB2(nearest);
ATH_MSG_DEBUG("Found distance of " << nearest << " from BCID " << bcid
<< " to the nearest unpaired B2 ");
// Add 127 to current BCID to check if inbound B1 is nearby FASER
- int offset_bcid = (bcid + 127) % 3564;
- nearest = findNearest(offset_bcid, bcid_mask, 1); // Inbound B1
- int nearest2 = findNearest(offset_bcid, bcid_mask, 3); // Inbound B1
- if (abs(nearest2) < abs(nearest)) nearest = nearest2;
+ if (m_lhcTool->getBeam1Bunches(ctx) == 0) {
+ ATH_MSG_INFO("No beam 1 bunches, can't set nearest");
+ nearest = -3564;
+ } else {
+ int offset_bcid = (bcid + 127) % 3564;
+ nearest = findNearest(offset_bcid, bcid_mask, 1); // Inbound B1
+ int nearest2 = findNearest(offset_bcid, bcid_mask, 3); // Inbound B1
+ if (abs(nearest2) < abs(nearest)) nearest = nearest2;
+ }
lhcDataHandle->set_distanceToInboundB1(nearest);
ATH_MSG_DEBUG("Found distance of " << nearest << " from BCID " << bcid
<< " to the nearest inbound B1 ");
unsigned int previous;
- previous = previousColliding(bcid, bcid_mask);
+ if (m_lhcTool->getCollidingBunches(ctx) == 0) {
+ ATH_MSG_INFO("No colliding bunches, can't set nearest");
+ previous = 9999;
+ } else {
+ previous = previousColliding(bcid, bcid_mask);
+ }
lhcDataHandle->set_distanceToPreviousColliding(previous);
ATH_MSG_DEBUG("Found distance of " << previous << " from BCID " << bcid
<< " to the previous colliding bunch ");
- previous = previousTrain(bcid, bcid_mask);
+ if (m_lhcTool->getCollidingBunches(ctx) == 0) {
+ ATH_MSG_INFO("No colliding bunches, can't set nearest");
+ previous = 9999;
+ } else {
+ previous = previousTrain(bcid, bcid_mask);
+ }
lhcDataHandle->set_distanceToTrainStart(previous);
ATH_MSG_DEBUG("Found distance of " << previous << " from BCID " << bcid
<< " to the start of the previous train ");
@@ -142,8 +177,20 @@ LHCDataAlg::findNearest(unsigned int bcid, const std::vector<unsigned char>& bci
// If we got to here, there was no match
// Does the BCID make sense?
- ATH_MSG_WARNING("Couldn't find distance from BCID " << bcid << " and pattern " << mask << "!");
- ATH_MSG_WARNING(bcid_mask);
+ ATH_MSG_WARNING("Couldn't find distance from BCID " << bcid << " and pattern " << int(mask) << "!");
+ int b1=0;
+ int b2=0;
+ int col = 0;
+ for (int i=0; i<3564; i++) {
+ if (mask & 0x01) b1++;
+ if (mask & 0x02) b2++;
+ if (mask & 0x03) col++;
+
+ ATH_MSG_WARNING("BCID " << i << " - " << int(bcid_mask[i]));
+ }
+
+ //ATH_MSG_WARNING(bcid_mask);
+ ATH_MSG_WARNING("Found B1: " << b1 << " B2: " << b2 << " Coll: " << col);
return -3565;
}
diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
index fab78641d7ae2a3842904122968a80170bca9a28..607c58d892d201731098f64829a82d556f1a302a 100644
--- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -1,11 +1,14 @@
atlas_subdir( NtupleDumper )
+find_package( nlohmann_json )
+
atlas_add_component(
NtupleDumper
src/NtupleDumperAlg.h
src/NtupleDumperAlg.cxx
src/component/NtupleDumper_entries.cxx
LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib
+ PRIVATE_LINK_LIBRARIES nlohmann_json::nlohmann_json
)
atlas_install_python_modules( python/*.py )
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index 27b30e25df16dd8baf1c4ce4b82b6a9e6a0f7b6f..1f4b990c2fb7bbdef7e109898e91bf5e0171d49d 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -44,8 +44,20 @@ parser.add_argument("--isMC", action='store_true',
parser.add_argument("--partial", action='store_true',
help="Allow partial input files")
+parser.add_argument("--trigFilt", action='store_true',
+ help="apply trigger event filter")
+parser.add_argument("--scintFilt", action='store_true',
+ help="apply scintillator event filter")
+parser.add_argument("--NoTrackFilt", action='store_true',
+ help="Don't apply track event filter (default: do)")
+
parser.add_argument("--unblind", action='store_true',
help="Don't apply signal blinding (default: do)")
+parser.add_argument("--onlyblind", action='store_true',
+ help="Only store events that were blinded (will override unblind arg)")
+
+parser.add_argument("--grl", default="",
+ help="Specify GRL to apply")
parser.add_argument("--fluka", action='store_true',
help="Add FLUKA weights to ntuple")
@@ -138,8 +150,7 @@ if filepath.is_dir():
print(f"Run = {runstr}")
print(f"First = {firstseg}")
print(f"Last = {lastseg}")
- print(f"Args = {args.tag}")
- print(f"Blind = {not args.unblind}")
+ print(f"Tag = {args.tag}")
# Find any tags
tagstr = firststem.replace(f"{firstfaser}-{firstshort}-{runstr}-{firstseg}", "")
@@ -199,6 +210,13 @@ if len(args.outfile) > 0:
print("Output file:")
print(outfile)
+print(f"Trigger Filter = {args.trigFilt}")
+print(f"Scintillator Filter = {args.scintFilt}")
+print(f"Track Filter = {not args.NoTrackFilt}")
+print(f"Blind = {not args.unblind}")
+print(f"OnlyBlinded = {args.onlyblind}")
+print(f"GRL = {args.grl}")
+
# OK, lets run the job here
from AthenaCommon.Logging import log, logging
from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
@@ -232,18 +250,25 @@ ConfigFlags.lock()
acc = MainServicesCfg(ConfigFlags)
acc.merge(PoolReadCfg(ConfigFlags))
+# Create kwargs for NtupleDumper
+grl_kwargs = {}
+if args.grl:
+ grl_kwargs['GoodRunsList'] = args.grl
+ grl_kwargs['ApplyGoodRunsList'] = True
+
+mc_kwargs = {}
+if args.genie:
+ mc_kwargs['UseGenieWeights'] = True
+if args.fluka:
+ mc_kwargs['UseFlukaWeights'] = True
+
# algorithm
from NtupleDumper.NtupleDumperConfig import NtupleDumperAlgCfg
if args.isMC:
- 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))
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, **mc_kwargs))
else:
- acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind)))
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = (not args.NoTrackFilt), DoTrigFilter = args.trigFilt, **grl_kwargs))
if not args.verbose:
from AthenaConfiguration.ComponentFactory import CompFactory
@@ -251,6 +276,10 @@ if not args.verbose:
AthenaEventLoopMgr.EventPrintoutInterval=1000
acc.addService(AthenaEventLoopMgr)
+else:
+ nd = acc.getEventAlgo("NtupleDumperAlg")
+ nd.OutputLevel = VERBOSE
+
# Hack to avoid problem with our use of MC databases when isMC = False
if not args.isMC:
replicaSvc = acc.getService("DBReplicaSvc")
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 3f476edb40480e5871a2df39a3cd2695ea894e4c..59a2be962c55923976f578c7e403cd2aae64a899 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -18,9 +18,6 @@
#include <TH1F.h>
#include <numeric>
-constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
-
-
NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
ISvcLocator *pSvcLocator)
: AthReentrantAlgorithm(name, pSvcLocator),
@@ -123,6 +120,19 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_spacePointContainerKey.initialize());
+ // Read GRL if requested
+ if (!m_goodRunsList.value().empty()) {
+ ATH_MSG_INFO("Opening GRL file " << m_goodRunsList);
+ ATH_MSG_INFO(m_applyGoodRunsList);
+ std::ifstream f(m_goodRunsList);
+ m_grl = nlohmann::json::parse(f);
+ f.close();
+ ATH_MSG_INFO("Read GRL with " << m_grl.size() << " records");
+ } else {
+ // Make sure this is empty
+ m_grl.clear();
+ }
+
if (m_useFlukaWeights)
{
m_baseEventCrossSection = (m_flukaCrossSection * kfemtoBarnsPerMilliBarn)/m_flukaCollisions;
@@ -143,6 +153,7 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("eventID", &m_event_number, "eventID/I");
m_tree->Branch("eventTime", &m_event_time, "eventTime/I");
m_tree->Branch("BCID", &m_bcid, "BCID/I");
+ m_tree->Branch("inGRL", &m_in_grl, "inGRL/I");
m_tree->Branch("fillNumber", &m_fillNumber, "fillNumber/I");
m_tree->Branch("betaStar", &m_betaStar, "betaStar/F");
@@ -167,10 +178,15 @@ StatusCode NtupleDumperAlg::initialize()
addWaveBranches("Preshower",2,12);
addWaveBranches("Calo",4,0);
+ m_tree->Branch("ScintHit", &m_scintHit);
+
addCalibratedBranches("Calo",4,0);
addBranch("Calo_total_nMIP", &m_calo_total_nMIP);
addBranch("Calo_total_E_dep", &m_calo_total_E_dep);
addBranch("Calo_total_E_EM", &m_calo_total_E_EM);
+ addBranch("Calo_total_raw_E_EM", &m_calo_total_raw_E_EM);
+ addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged);
+ addBranch("Calo_total_raw_E_EM_fudged", &m_calo_total_raw_E_EM_fudged);
addCalibratedBranches("Preshower",2,12);
addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP);
@@ -261,6 +277,11 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("Track_ThetaX_atCalo", &m_thetaxCalo);
m_tree->Branch("Track_ThetaY_atCalo", &m_thetayCalo);
+ m_tree->Branch("Track_x_atMaxRadius", &m_x_atMaxRadius);
+ m_tree->Branch("Track_y_atMaxRadius", &m_y_atMaxRadius);
+ m_tree->Branch("Track_z_atMaxRadius", &m_z_atMaxRadius);
+ m_tree->Branch("Track_r_atMaxRadius", &m_r_atMaxRadius);
+
//TRUTH
m_tree->Branch("t_pdg", &m_t_pdg);
m_tree->Branch("t_barcode", &m_t_barcode);
@@ -385,10 +406,10 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14]));
- m_MIP_sim_Edep_calo = 0.0585; // MIP deposits 0.0585 GeV of energy in calo
- m_MIP_sim_Edep_preshower = 0.004894; // MIP deposits 0.004894 GeV of energy in a preshower layer
-
- if (m_doBlinding) {
+ if (m_onlyBlinded){
+ ATH_MSG_INFO("Only events that would be blinded are saved in ntuple");
+ m_doBlinding = false;
+ } else if (m_doBlinding) {
ATH_MSG_INFO("Blinding will be enforced for real data.");
} else {
ATH_MSG_INFO("Blinding will NOT be enforced for real data.");
@@ -410,85 +431,162 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
isMC = true;
}
- // if real data, store charge in histograms from random events and only fill ntuple from coincidence events
- if (!isMC) {
- SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx);
- m_tap=triggerData->tap();
- // unpack trigger word: 1=calo, 2=veotnu|veto1|preshower, 4=TimingLayer, 8=(VetoNu|Veto2)&Preshower, 16=random, 32=LED
- bool trig_random = false;
- if ( m_tap == 16 ) trig_random = true;
+ // EventInfo data
+ m_run_number = ctx.eventID().run_number();
+ m_event_number = ctx.eventID().event_number();
+ m_event_time = ctx.eventID().time_stamp();
+ m_bcid = ctx.eventID().bunch_crossing_id();
+
+ // For real data, find if data is in GoodRunsList
+ if (!isMC && !m_grl.empty()) {
+ ATH_MSG_DEBUG("Checking GRL for run " << m_run_number << " event " << m_event_number);
+
+ // JSON keys are always strings
+ std::string run_string = std::to_string(m_run_number);
+ std::string reason("unknown");
+
+ if (m_grl.contains(run_string)) {
+ // m_in_grl set to 0 by clearTree
+ auto jrun = m_grl[run_string];
+
+ // Make sure this is in the stable beams range
+ if (jrun.contains("stable_list")) {
+ for ( auto& slist : jrun["stable_list"]) {
+ // Check if event falls in this range
+ if (slist["start_utime"] > m_event_time) continue;
+ if (slist["stop_utime"] <= m_event_time) continue;
+ m_in_grl = 1;
+ break;
+ }
+
+ // Check exclude list (not all runs will have this)
+ if (jrun.contains("excluded_list")) {
+ ATH_MSG_DEBUG("Checking excluded list");
+ for ( auto& blist : jrun["excluded_list"]) {
+ ATH_MSG_DEBUG("Checking excluded list start " << blist["start_utime"] << " stop " << blist["stop_utime"] << " event " << m_event_time);
+
+ // Check if event falls in this range
+ if (blist["start_utime"] > m_event_time) continue;
+ if (blist["stop_utime"] <= m_event_time) continue;
+ if (blist.contains("reason")) {
+ reason = blist["reason"];
+ }
+ m_in_grl = 0;
+ break;
+ }
+ }
- bool trig_coincidence_preshower_and_vetoes = false;
- if ( (m_tap&8) != 0 ) trig_coincidence_preshower_and_vetoes = true;
+ if (m_in_grl == 0) {
+ ATH_MSG_DEBUG("Run " << run_string << " event " << m_event_number << " time " << m_event_time << " excluded for " << reason);
+ }
- bool trig_coincidence_timing_and_vetoesORpreshower = false;
- if ( ((m_tap&4)!=0) && ((m_tap&2)!=0) ) trig_coincidence_timing_and_vetoesORpreshower = true;
+ } else {
+ ATH_MSG_WARNING("Run " << run_string << " has no stable_list!");
+ }
+ } else {
+ ATH_MSG_WARNING("Run " << run_string << " not found in GRL!");
+ }
+ }
- bool trig_coincidence_timing_and_calo = false;
- if ( ((m_tap&4)!=0) && ((m_tap&1)!=0) ) trig_coincidence_timing_and_calo = true;
+ // Skip this event?
+ if (m_applyGoodRunsList && (m_in_grl == 0)) {
+ ATH_MSG_DEBUG("Skipping run " << m_run_number
+ << " event " << m_event_number << " - outside GRL" );
+ m_eventsFailedGRL += 1;
+ return StatusCode::SUCCESS;
+ }
- bool trig_coincidence_vetoesORpreshower_and_calo = false;
- if ( ((m_tap&2)!=0) && ((m_tap&1)!=0) ) trig_coincidence_vetoesORpreshower_and_calo = true;
+ // if real data, correct for diigitzer timing jitter with clock phase
+ if (!isMC) {
+ // correct waveform time with clock phase
+ // needs to be done before processing waveform values
+ SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_ClockWaveformContainer, ctx);
+ ATH_CHECK(clockHandle.isValid());
- bool trig_calo = (m_tap & 1);
+ if (clockHandle->phase() < -2.0) { // wrap around clock pahse so -pi goes to pi
+ m_clock_phase = ((clockHandle->phase() + 2*3.14159) / 3.14159) * 12.5;
+ } else {
+ m_clock_phase = (clockHandle->phase() / 3.14159) * 12.5;
+ }
+ }
- // for random trigger, store charge of scintillators in histograms
- if (trig_random) {
- // Read in Waveform containers
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
- ATH_CHECK(vetoNuContainer.isValid());
+ // process all waveform data for all scintillator and calorimeter channels
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
+ ATH_CHECK(vetoNuContainer.isValid());
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
- ATH_CHECK(vetoContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
+ ATH_CHECK(vetoContainer.isValid());
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ ATH_CHECK(triggerContainer.isValid());
- SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
- ATH_CHECK(triggerContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
+ ATH_CHECK(preshowerContainer.isValid());
- SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
- ATH_CHECK(preshowerContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
+ ATH_CHECK(ecalContainer.isValid());
- SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
- ATH_CHECK(ecalContainer.isValid());
+ FillWaveBranches(*vetoNuContainer);
+ FillWaveBranches(*vetoContainer);
+ FillWaveBranches(*triggerContainer);
+ FillWaveBranches(*preshowerContainer);
+ FillWaveBranches(*ecalContainer);
+ // if real data, store charge in histograms from random events
+ if (!isMC) {
+ SG::ReadHandle<xAOD::FaserTriggerData> triggerData(m_FaserTriggerData, ctx);
+ m_tap=triggerData->tap();
+ // for random trigger, store charge of scintillators in histograms
+ if ((m_tap&16) != 0) {
// Fill histograms
- if (vetoNuContainer.isValid()) {
- for (auto hit : *vetoNuContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (vetoContainer.isValid()) {
- for (auto hit : *vetoContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (triggerContainer.isValid()) {
- for (auto hit : *triggerContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (preshowerContainer.isValid()) {
- for (auto hit : *preshowerContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
- }
- if (ecalContainer.isValid()) {
- for (auto hit : *ecalContainer) {
- int ch=hit->channel();
- m_HistRandomCharge[ch]->Fill(hit->raw_integral()/50.0);
- }
+ for (int chan = 0; chan<15; chan++) {
+ m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]);
}
+ return StatusCode::SUCCESS; // finished with this randomly triggered event
+ }
+
+ if (m_doTrigFilter) {
+
+ bool trig_coincidence_preshower_and_vetoes = ( (m_tap&8) != 0 );
+ bool trig_coincidence_timing_and_vetoesORpreshower = ( ((m_tap&4)!=0) && ((m_tap&2)!=0) );
+ bool trig_calo = (m_tap & 1);
- return StatusCode::SUCCESS; // finished with this randomly triiggered event
+ if ( !(trig_coincidence_preshower_and_vetoes || trig_coincidence_timing_and_vetoesORpreshower || trig_calo) ) {
+ // don't process events that fail to activate coincidence triggers
+ ATH_MSG_DEBUG("event did not pass trigger filter");
+ return StatusCode::SUCCESS;
+ }
+ }
- // 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;
+ if (m_doScintFilter) {
+ // filter events, but use waveform peak cuts instead of triggers, as triggers could miss signals slightly out of time
+ // digitizer channels described here: https://twiki.cern.ch/twiki/bin/viewauth/FASER/CaloScintillator
+ bool calo_trig = (m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0);
+ bool vetoNu_trig = (m_wave_raw_peak[4] > 25.0 && m_wave_raw_peak[5] > 25.0);
+ bool vetoSt2_trig = (m_wave_raw_peak[6] > 25.0 && m_wave_raw_peak[7] > 25.0);
+ bool timingScint_trig = ((m_wave_raw_peak[8] > 25.0 && m_wave_raw_peak[9] > 25.0) || (m_wave_raw_peak[10] > 25.0 && m_wave_raw_peak[11] > 25.0));
+ bool preshower_trig = (m_wave_raw_peak[12] > 3.0 && m_wave_raw_peak[13] > 3.0);
+ bool vetoSt1_trig = (m_wave_raw_peak[14] > 25.0);
+
+ bool veto_OR_trig = (vetoNu_trig || vetoSt1_trig || vetoSt2_trig);
+
+ bool passes_ScintFilter = false;
+ if (calo_trig) {
+ passes_ScintFilter = true;
+ } else if (veto_OR_trig && timingScint_trig) {
+ passes_ScintFilter = true;
+ } else if (veto_OR_trig && preshower_trig) {
+ passes_ScintFilter = true;
+ } else if (timingScint_trig && preshower_trig) {
+ passes_ScintFilter = true;
+ }
+
+ if (!passes_ScintFilter) {
+ ATH_MSG_DEBUG("event did not pass scint filter");
+ return StatusCode::SUCCESS; // only store events that pass filter
+ }
}
+
// store trigger data in ntuple variables
m_tbp=triggerData->tbp();
m_tap=triggerData->tap();
@@ -530,144 +628,149 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("LHC data distanceToTrainStart = " << lhcData->distanceToTrainStart() );
ATH_MSG_DEBUG("LHC data distanceToPreviousColliding = " << lhcData->distanceToPreviousColliding() );
- // correct waveform time with clock phase
- SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_ClockWaveformContainer, ctx);
- ATH_CHECK(clockHandle.isValid());
-
- if (clockHandle->phase() < -2.0) { // wrap around clock pahse so -pi goes to pi
- m_clock_phase = ((clockHandle->phase() + 2*3.14159) / 3.14159) * 12.5;
- } else {
- m_clock_phase = (clockHandle->phase() / 3.14159) * 12.5;
- }
-
} // done with processing only on real data
- m_run_number = ctx.eventID().run_number();
- m_event_number = ctx.eventID().event_number();
- m_event_time = ctx.eventID().time_stamp();
- m_bcid = ctx.eventID().bunch_crossing_id();
+ // fill scintHit word with bits that reflect if a scintillator was hit (1 = vetoNu0, 2 = vetoNu1, 4 = vetoSt1_1, 8 vetoSt2_0, 16 = vetoSt2_1, 32 = Timing scint, 64 = preshower0, 128 = preshower1)
+ if (m_wave_raw_charge[4] > 40.0) {
+ m_scintHit = m_scintHit | 1;
+ }
+ if (m_wave_raw_charge[5] > 40.0) {
+ m_scintHit = m_scintHit | 2;
+ }
+ if (m_wave_raw_charge[14] > 40.0) {
+ m_scintHit = m_scintHit | 4;
+ }
+ if (m_wave_raw_charge[6] > 40.0) {
+ m_scintHit = m_scintHit | 8;
+ }
+ if (m_wave_raw_charge[7] > 40.0) {
+ m_scintHit = m_scintHit | 16;
+ }
+ if (m_wave_raw_charge[8]+m_wave_raw_charge[9]+m_wave_raw_charge[10]+m_wave_raw_charge[11] > 25.0) {
+ m_scintHit = m_scintHit | 32;
+ }
+ if (m_wave_raw_charge[12] > 2.5) {
+ m_scintHit = m_scintHit | 64;
+ }
+ if (m_wave_raw_charge[13] > 2.5) {
+ m_scintHit = m_scintHit | 128;
+ }
if (isMC) { // if simulation find MC cross section and primary lepton
// Work out effective cross section for MC
- if (m_useFlukaWeights)
- {
- double flukaWeight = truthEventContainer->at(0)->weights()[0];
- ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
- m_crossSection = m_baseEventCrossSection * flukaWeight;
- }
- else if (m_useGenieWeights)
- {
- m_crossSection = m_baseEventCrossSection;
- }
- else
- {
+ if (m_useFlukaWeights) {
+ double flukaWeight = truthEventContainer->at(0)->weights()[0];
+ ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
+ m_crossSection = m_baseEventCrossSection * flukaWeight;
+ } else if (m_useGenieWeights) {
+ m_crossSection = m_baseEventCrossSection;
+ } else {
//ATH_MSG_WARNING("Monte carlo event with no weighting scheme specified. Setting crossSection (weight) to " << m_baseEventCrossSection << " fb.");
- m_crossSection = m_baseEventCrossSection;
+ m_crossSection = m_baseEventCrossSection;
}
- // Find the M d0 and d1 truth information
+ // Find truth particle information
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
- if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
- {
-
+ if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) {
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());
- m_truthd0_px.push_back(particle->p4().X());
- m_truthd0_py.push_back(particle->p4().Y());
- m_truthd0_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthd0_x.push_back(particle->prodVtx()->x());
- m_truthd0_y.push_back(particle->prodVtx()->y());
- m_truthd0_z.push_back(particle->prodVtx()->z());
- } else {
- m_truthd0_x.push_back(NaN);
- m_truthd0_y.push_back(NaN);
- m_truthd0_z.push_back(NaN);
- }
- }
- if ( particle->pdgId() == -11) // daughter particle (electron)
- {
- m_truthd1_P.push_back(particle->p4().P());
- m_truthd1_px.push_back(particle->p4().X());
- m_truthd1_py.push_back(particle->p4().Y());
- m_truthd1_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthd1_x.push_back(particle->prodVtx()->x());
- m_truthd1_y.push_back(particle->prodVtx()->y());
- m_truthd1_z.push_back(particle->prodVtx()->z());
- } else {
- m_truthd1_x.push_back(NaN);
- m_truthd1_y.push_back(NaN);
- m_truthd1_z.push_back(NaN);
- }
- }
- }
+ for (auto particle : *truthParticleContainer) { // loop over first 10 truth particles (for non A' samples)
+ if (ipart++ > 9) break;
+
+ 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(999999);
+ m_truth_prod_y.push_back(999999);
+ m_truth_prod_z.push_back(999999);
+ }
+
+ 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(999999);
+ m_truth_dec_y.push_back(999999);
+ m_truth_dec_z.push_back(999999);
+ }
+
+ // Find the M d0 and d1 truth information for dark photon
+ if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
+ auto positions = m_fiducialParticleTool->getTruthPositions(particle->barcode());
+ 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(999999);
+ m_truthM_y.push_back(999999);
+ m_truthM_z.push_back(999999);
+ }
+ }
+
+ if ( particle->pdgId() == 11) { // daughter particle (electron)
+ m_truthd0_P.push_back(particle->p4().P());
+ m_truthd0_px.push_back(particle->p4().X());
+ m_truthd0_py.push_back(particle->p4().Y());
+ m_truthd0_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasProdVtx()) {
+ m_truthd0_x.push_back(particle->prodVtx()->x());
+ m_truthd0_y.push_back(particle->prodVtx()->y());
+ m_truthd0_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd0_x.push_back(999999);
+ m_truthd0_y.push_back(999999);
+ m_truthd0_z.push_back(999999);
+ }
+ for (int station = 1; station < 4; ++station) {
+ m_truthd0_x.push_back(positions[station].x());
+ m_truthd0_y.push_back(positions[station].y());
+ m_truthd0_z.push_back(positions[station].z());
+ }
+ }
+ if ( particle->pdgId() == -11) { // daughter particle (positron)
+ m_truthd1_P.push_back(particle->p4().P());
+ m_truthd1_px.push_back(particle->p4().X());
+ m_truthd1_py.push_back(particle->p4().Y());
+ m_truthd1_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasProdVtx()) {
+ m_truthd1_x.push_back(particle->prodVtx()->x());
+ m_truthd1_y.push_back(particle->prodVtx()->y());
+ m_truthd1_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd1_x.push_back(999999);
+ m_truthd1_y.push_back(999999);
+ m_truthd1_z.push_back(999999);
+ }
+ for (int station = 1; station < 4; ++station) {
+ m_truthd1_x.push_back(positions[station].x());
+ m_truthd1_y.push_back(positions[station].y());
+ m_truthd1_z.push_back(positions[station].z());
+ }
+ }
+ }
}
}
}
// load in calibrated calo container
- if (m_doCalib) {
SG::ReadHandle<xAOD::CalorimeterHitContainer> ecalCalibratedContainer { m_ecalCalibratedContainer, ctx };
ATH_CHECK(ecalCalibratedContainer.isValid());
for (auto hit : *ecalCalibratedContainer) {
@@ -679,6 +782,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_calo_total_nMIP += hit->Nmip();
m_calo_total_E_dep += hit->E_dep();
m_calo_total_E_EM += hit->E_EM();
+ m_calo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio());
ATH_MSG_DEBUG("Calibrated calo: 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());
@@ -694,6 +798,14 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
//}
}
+ // add a fudged energy variable that corrects the MC to agree with the testbeam results
+ m_calo_total_E_EM_fudged = m_calo_total_E_EM;
+ m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM;
+ if (isMC) {
+ m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
+ m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM_fudged * m_caloMC_FudgeFactor;
+ }
+
// load in calibrated preshower container
SG::ReadHandle<xAOD::CalorimeterHitContainer> preshowerCalibratedContainer { m_preshowerCalibratedContainer, ctx };
ATH_CHECK(preshowerCalibratedContainer.isValid());
@@ -707,40 +819,24 @@ 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 };
- ATH_CHECK(vetoNuContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
- ATH_CHECK(vetoContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
- ATH_CHECK(triggerContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
- ATH_CHECK(preshowerContainer.isValid());
-
- SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
- ATH_CHECK(ecalContainer.isValid());
-
- FillWaveBranches(*vetoNuContainer);
- FillWaveBranches(*vetoContainer);
- FillWaveBranches(*triggerContainer);
- FillWaveBranches(*preshowerContainer);
- 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
// Only blind colliding BCIDs (+/- 1)
- if ( (!isMC) && m_doBlinding && abs(m_distanceToCollidingBCID) <= 1) {
+ bool blinded_event = false;
+ if ((!isMC) && 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;
+ if (m_wave_raw_charge[4] < 40.0 and m_wave_raw_charge[5] < 40.0 and m_wave_raw_charge[6] < 40.0 and m_wave_raw_charge[7] < 40.0 and m_wave_raw_charge[14] < 40.0) { // channels 4 and 5 are vetoNu, channels 6,7, and 14 are veto
+ blinded_event = true;
}
}
}
+ if (blinded_event && m_doBlinding) {
+ return StatusCode::SUCCESS;
+ } else if (!blinded_event && m_onlyBlinded) {
+ 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;
@@ -753,73 +849,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext();
auto gctx = faserGeometryContext.context();
- // loop over clusters and store how many clusters are in each tracking station
- SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
- ATH_CHECK(clusterContainer.isValid());
-
- for (auto collection : *clusterContainer)
- {
- Identifier id = collection->identify();
- int station = m_sctHelper->station(id);
- int clusters = (int) collection->size();
- switch (station)
- {
- case 0:
- m_station0Clusters += clusters;
- // following lines commented out depict how to access cluster position
- //for (auto cluster : *collection) {
- // if (cluster == nullptr) continue;
- // auto pos = cluster->globalPosition();
- // m_station0ClusterX.push_back(pos.x());
- //}
- break;
- case 1:
- m_station1Clusters += clusters;
- break;
- case 2:
- m_station2Clusters += clusters;
- break;
- case 3:
- m_station3Clusters += clusters;
- break;
- default:
- ATH_MSG_FATAL("Unknown tracker station number " << station);
- break;
- }
- }
-
- // loop over spacepoints and store each space point position
- SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
- ATH_CHECK(spacePointContainer.isValid());
- for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
- m_nspacepoints += spacePointCollection->size();
- for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
- auto pos = spacePoint->globalPosition();
- m_spacepointX.push_back(pos.x());
- m_spacepointY.push_back(pos.y());
- m_spacepointZ.push_back(pos.z());
- }
- }
-
- // loop over track segments and store position, momentum, chi2, and nDOF for each segment
- SG::ReadHandle<TrackCollection> trackSegmentCollection {m_trackSegmentCollection, ctx};
- ATH_CHECK(trackSegmentCollection.isValid());
- for (const Trk::Track* trackSeg : *trackSegmentCollection) {
- if (trackSeg == nullptr) continue;
- m_ntracksegs += 1;
- m_trackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
- m_trackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
- auto SegParameters = trackSeg->trackParameters()->front();
- const Amg::Vector3D SegPosition = SegParameters->position();
- const Amg::Vector3D SegMomentum = SegParameters->momentum();
- m_trackseg_x.push_back(SegPosition.x());
- m_trackseg_y.push_back(SegPosition.y());
- m_trackseg_z.push_back(SegPosition.z());
- m_trackseg_px.push_back(SegMomentum.x());
- m_trackseg_py.push_back(SegMomentum.y());
- m_trackseg_pz.push_back(SegMomentum.z());
- }
-
// Write out all truth particle barcodes that have a momentum larger than MinMomentum (default is 50 GeV)
std::map<int, size_t> truthParticleCount {};
if (isMC) {
@@ -836,7 +865,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// extrapolate track to all scintillator positions and store extrapolated position and angle
SG::ReadHandle<TrackCollection> trackCollection;
if (m_useIFT) {
- SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that excludes IFT
+ SG::ReadHandle<TrackCollection> tc {m_trackCollection, ctx}; // use track collection that includes IFT
trackCollection = tc;
} else {
SG::ReadHandle<TrackCollection> tc {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
@@ -896,6 +925,20 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_pzdown.push_back(downstreamParameters->momentum().z());
m_pdown.push_back(sqrt( pow(downstreamParameters->momentum().x(),2) + pow(downstreamParameters->momentum().y(),2) + pow(downstreamParameters->momentum().z(),2) ));
+ // find and store the position of the measurement on track with the largest radius
+ double maxRadius = 0;
+ Amg::Vector3D positionAtMaxRadius {};
+ for (const Trk::TrackParameters* trackParameters : *track->trackParameters()) {
+ if (radius(trackParameters->position()) > maxRadius) {
+ maxRadius = radius(trackParameters->position());
+ positionAtMaxRadius = trackParameters->position();
+ }
+ }
+ m_r_atMaxRadius.push_back(maxRadius);
+ m_x_atMaxRadius.push_back(positionAtMaxRadius.x());
+ m_y_atMaxRadius.push_back(positionAtMaxRadius.y());
+ m_z_atMaxRadius.push_back(positionAtMaxRadius.z());
+
if (isMC) { // if simulation, store track truth info as well
auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
if (truthParticle != nullptr) {
@@ -917,18 +960,18 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
m_t_prodVtx_z.push_back(truthParticle->prodVtx()->z());
} else {
- m_t_prodVtx_x.push_back(NaN);
- m_t_prodVtx_y.push_back(NaN);
- m_t_prodVtx_z.push_back(NaN);
+ m_t_prodVtx_x.push_back(999999);
+ m_t_prodVtx_y.push_back(999999);
+ m_t_prodVtx_z.push_back(999999);
}
if (truthParticle->hasDecayVtx()) {
m_t_decayVtx_x.push_back(truthParticle->decayVtx()->x());
m_t_decayVtx_y.push_back(truthParticle->decayVtx()->y());
m_t_decayVtx_z.push_back(truthParticle->decayVtx()->z());
} else {
- m_t_decayVtx_x.push_back(NaN);
- m_t_decayVtx_y.push_back(NaN);
- m_t_decayVtx_z.push_back(NaN);
+ m_t_decayVtx_x.push_back(999999);
+ m_t_decayVtx_y.push_back(999999);
+ m_t_decayVtx_z.push_back(999999);
}
m_t_px.push_back(truthParticle->px());
m_t_py.push_back(truthParticle->py());
@@ -940,167 +983,175 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_t_eta.push_back(truthParticle->p4().Eta());
} else {
ATH_MSG_WARNING("Can not find truthParticle.");
- setNaN();
+ clearTrackTruth();
}
} else {
- setNaN();
+ clearTrackTruth();
}
- // fill extrapolation vectors with NaN, will get set to real number if the track extrapolation succeeds
- m_xVetoNu.push_back(NaN);
- m_yVetoNu.push_back(NaN);
- m_thetaxVetoNu.push_back(NaN);
- m_thetayVetoNu.push_back(NaN);
- m_xVetoStation1.push_back(NaN);
- m_yVetoStation1.push_back(NaN);
- m_thetaxVetoStation1.push_back(NaN);
- m_thetayVetoStation1.push_back(NaN);
- m_xVetoStation2.push_back(NaN);
- m_yVetoStation2.push_back(NaN);
- m_thetaxVetoStation2.push_back(NaN);
- m_thetayVetoStation2.push_back(NaN);
- m_xTrig.push_back(NaN);
- m_yTrig.push_back(NaN);
- m_thetaxTrig.push_back(NaN);
- m_thetayTrig.push_back(NaN);
- m_xPreshower1.push_back(NaN);
- m_yPreshower1.push_back(NaN);
- m_thetaxPreshower1.push_back(NaN);
- m_thetayPreshower1.push_back(NaN);
- m_xPreshower2.push_back(NaN);
- m_yPreshower2.push_back(NaN);
- m_thetaxPreshower2.push_back(NaN);
- m_thetayPreshower2.push_back(NaN);
- m_xCalo.push_back(NaN);
- m_yCalo.push_back(NaN);
- m_thetaxCalo.push_back(NaN);
- m_thetayCalo.push_back(NaN);
-
- // extrapolate track from first station
- if (stationMap.count(1) > 0) { // extrapolation crashes if the track parameters are is too far away to extrapolate
- Amg::Vector3D position = upstreamParameters->position();
- Amg::Vector3D momentum = upstreamParameters->momentum();
- Acts::BoundVector params = Acts::BoundVector::Zero();
- params[Acts::eBoundLoc0] = -position.y();
- params[Acts::eBoundLoc1] = position.x();
- params[Acts::eBoundPhi] = momentum.phi();
- params[Acts::eBoundTheta] = momentum.theta();
- params[Acts::eBoundQOverP] = upstreamParameters->charge() / momentum.mag();
- params[Acts::eBoundTime] = 0;
- auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParameters(std::move(startSurface), params, upstreamParameters->charge());
-
- auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_VetoNu, Acts::backward);
- if (targetParameters_VetoNu != nullptr) {
- auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
- auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
- m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
- m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
- m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
- m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
- } else {
- ATH_MSG_INFO("vetoNu null targetParameters");
- }
-
- auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto1, Acts::backward);
- if (targetParameters_Veto1 != nullptr) {
- auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
- auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
- m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
- m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
- m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
- m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
- } else {
- ATH_MSG_INFO("veto1 null targetParameters");
- }
+ // fill extrapolation vectors with dummy values, will get set to real number if the track extrapolation succeeds
+ m_xVetoNu.push_back(999999);
+ m_yVetoNu.push_back(999999);
+ m_thetaxVetoNu.push_back(999999);
+ m_thetayVetoNu.push_back(999999);
+ m_xVetoStation1.push_back(999999);
+ m_yVetoStation1.push_back(999999);
+ m_thetaxVetoStation1.push_back(999999);
+ m_thetayVetoStation1.push_back(999999);
+ m_xVetoStation2.push_back(999999);
+ m_yVetoStation2.push_back(999999);
+ m_thetaxVetoStation2.push_back(999999);
+ m_thetayVetoStation2.push_back(999999);
+ m_xTrig.push_back(999999);
+ m_yTrig.push_back(999999);
+ m_thetaxTrig.push_back(999999);
+ m_thetayTrig.push_back(999999);
+ m_xPreshower1.push_back(999999);
+ m_yPreshower1.push_back(999999);
+ m_thetaxPreshower1.push_back(999999);
+ m_thetayPreshower1.push_back(999999);
+ m_xPreshower2.push_back(999999);
+ m_yPreshower2.push_back(999999);
+ m_thetaxPreshower2.push_back(999999);
+ m_thetayPreshower2.push_back(999999);
+ m_xCalo.push_back(999999);
+ m_yCalo.push_back(999999);
+ m_thetaxCalo.push_back(999999);
+ m_thetayCalo.push_back(999999);
+
+ // Define paramters for track extrapolation from most upstream measurement
+ Amg::Vector3D position_up = upstreamParameters->position();
+ Amg::Vector3D momentum_up = upstreamParameters->momentum();
+ Acts::BoundVector params_up = Acts::BoundVector::Zero();
+ params_up[Acts::eBoundLoc0] = -position_up.y();
+ params_up[Acts::eBoundLoc1] = position_up.x();
+ params_up[Acts::eBoundPhi] = momentum_up.phi();
+ params_up[Acts::eBoundTheta] = momentum_up.theta();
+ params_up[Acts::eBoundQOverP] = upstreamParameters->charge() / momentum_up.mag();
+ params_up[Acts::eBoundTime] = 0;
+ auto startSurface_up = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_up.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_up(std::move(startSurface_up), params_up, upstreamParameters->charge());
+
+ // Define paramters for track extrapolation from most downstream measurement
+ Amg::Vector3D position_down = downstreamParameters->position();
+ Amg::Vector3D momentum_down = downstreamParameters->momentum();
+ Acts::BoundVector params_down = Acts::BoundVector::Zero();
+ params_down[Acts::eBoundLoc0] = -position_down.y();
+ params_down[Acts::eBoundLoc1] = position_down.x();
+ params_down[Acts::eBoundPhi] = momentum_down.phi();
+ params_down[Acts::eBoundTheta] = momentum_down.theta();
+ params_down[Acts::eBoundQOverP] = downstreamParameters->charge() / momentum_down.mag();
+ params_down[Acts::eBoundTime] = 0;
+ auto startSurface_down = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_down.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_down(std::move(startSurface_down), params_down, downstreamParameters->charge());
+
+ // Extrapolate track to scintillators
+ auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_VetoNu, Acts::backward);
+ if (targetParameters_VetoNu != nullptr) {
+ auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
+ auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
+ m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
+ m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
+ m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
+ m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
+ } else {
+ ATH_MSG_INFO("vetoNu null targetParameters");
+ }
- auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Veto2, Acts::backward);
- if (targetParameters_Veto2 != nullptr) {
- auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
- auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
- m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
- m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
- m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
- m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
- } else {
- ATH_MSG_INFO("veto2 null targetParameters");
- }
+ auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto1, Acts::backward);
+ if (targetParameters_Veto1 != nullptr) {
+ auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
+ auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
+ m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
+ m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
+ m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
+ m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
+ } else {
+ ATH_MSG_INFO("veto1 null targetParameters");
+ }
- auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters, *targetSurface_Trig, Acts::backward); // must extrapolate backsards to trig plane if track starts in station 1
- if (targetParameters_Trig != nullptr) {
- auto targetPosition_Trig = targetParameters_Trig->position(gctx);
- auto targetMomentum_Trig = targetParameters_Trig->momentum();
- m_xTrig[m_longTracks] = targetPosition_Trig.x();
- m_yTrig[m_longTracks] = targetPosition_Trig.y();
- m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
- m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
- } else {
- ATH_MSG_INFO("Trig null targetParameters");
- }
+ auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto2, Acts::backward);
+ if (targetParameters_Veto2 != nullptr) {
+ auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
+ auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
+ m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
+ m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
+ m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
+ m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
+ } else {
+ ATH_MSG_INFO("veto2 null targetParameters");
+ }
+ auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Trig, Acts::backward); // must extrapolate backsards to trig plane if track starts in station 1
+ if (targetParameters_Trig != nullptr) {
+ auto targetPosition_Trig = targetParameters_Trig->position(gctx);
+ auto targetMomentum_Trig = targetParameters_Trig->momentum();
+ m_xTrig[m_longTracks] = targetPosition_Trig.x();
+ m_yTrig[m_longTracks] = targetPosition_Trig.y();
+ m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
+ m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
+ } else {
+ ATH_MSG_INFO("Trig null targetParameters");
}
- // extrapolate track from tracking station 3
- if (stationMap.count(3) > 0) { // extrapolation crashes if the track does not end in the Station 3, as it is too far away to extrapolate
- Amg::Vector3D positionDown = downstreamParameters->position();
- Amg::Vector3D momentumDown = downstreamParameters->momentum();
- Acts::BoundVector paramsDown = Acts::BoundVector::Zero();
- paramsDown[Acts::eBoundLoc0] = -positionDown.y();
- paramsDown[Acts::eBoundLoc1] = positionDown.x();
- paramsDown[Acts::eBoundPhi] = momentumDown.phi();
- paramsDown[Acts::eBoundTheta] = momentumDown.theta();
- paramsDown[Acts::eBoundQOverP] = downstreamParameters->charge() / momentumDown.mag();
- paramsDown[Acts::eBoundTime] = 0;
- auto startSurfaceDown = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, positionDown.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParametersDown(std::move(startSurfaceDown), paramsDown, downstreamParameters->charge());
-
- auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower1, Acts::forward);
- if (targetParameters_Preshower1 != nullptr) {
- auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
- auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
- m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
- m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
- m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
- m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
- } else {
- ATH_MSG_INFO("Preshower1 null targetParameters");
- }
+ auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower1, Acts::forward);
+ if (targetParameters_Preshower1 != nullptr) {
+ auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
+ auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
+ m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
+ m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
+ m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
+ m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
+ } else {
+ ATH_MSG_INFO("Preshower1 null targetParameters");
+ }
- auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Preshower2, Acts::forward);
- if (targetParameters_Preshower2 != nullptr) {
- auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
- auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
- m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
- m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
- m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
- m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
- } else {
- ATH_MSG_INFO("Preshower2 null targetParameters");
- }
+ auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower2, Acts::forward);
+ if (targetParameters_Preshower2 != nullptr) {
+ auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
+ auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
+ m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
+ m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
+ m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
+ m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
+ } else {
+ ATH_MSG_INFO("Preshower2 null targetParameters");
+ }
- auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
- std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParametersDown, *targetSurface_Calo, Acts::forward);
- if (targetParameters_Calo != nullptr) {
- auto targetPosition_Calo = targetParameters_Calo->position(gctx);
- auto targetMomentum_Calo = targetParameters_Calo->momentum();
- m_xCalo[m_longTracks] = targetPosition_Calo.x();
- m_yCalo[m_longTracks] = targetPosition_Calo.y();
- m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
- m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
- } else {
- ATH_MSG_INFO("Calo null targetParameters");
- }
+ auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Calo, Acts::forward);
+ if (targetParameters_Calo != nullptr) {
+ auto targetPosition_Calo = targetParameters_Calo->position(gctx);
+ auto targetMomentum_Calo = targetParameters_Calo->momentum();
+ m_xCalo[m_longTracks] = targetPosition_Calo.x();
+ m_yCalo[m_longTracks] = targetPosition_Calo.y();
+ m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
+ m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
+ } else {
+ ATH_MSG_INFO("Calo null targetParameters");
}
m_longTracks++;
}
+ if (!isMC) {
+ if (m_doTrackFilter) { // filter events: colliding bunches have at least one long track, non-colliding bunches have at least one long track or one calo module with raw_peak > 3 mV
+ if (m_longTracks == 0 && abs(m_distanceToCollidingBCID) <= 1) {
+ return StatusCode::SUCCESS;
+ } else if (abs(m_distanceToCollidingBCID) > 1) {
+ if (!(m_longTracks > 0 || m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0 )) {
+ return StatusCode::SUCCESS;
+ }
+ }
+ }
+ }
+
if (isMC) {
for (auto &tp : truthParticleCount) {
m_truthParticleBarcode.push_back(tp.first);
@@ -1109,40 +1160,99 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
}
+ // loop over clusters and store how many clusters are in each tracking station
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+ ATH_CHECK(clusterContainer.isValid());
+
+ for (auto collection : *clusterContainer)
+ {
+ Identifier id = collection->identify();
+ int station = m_sctHelper->station(id);
+ int clusters = (int) collection->size();
+ switch (station)
+ {
+ case 0:
+ m_station0Clusters += clusters;
+ // following lines commented out depict how to access cluster position
+ //for (auto cluster : *collection) {
+ // if (cluster == nullptr) continue;
+ // auto pos = cluster->globalPosition();
+ // m_station0ClusterX.push_back(pos.x());
+ //}
+ break;
+ case 1:
+ m_station1Clusters += clusters;
+ break;
+ case 2:
+ m_station2Clusters += clusters;
+ break;
+ case 3:
+ m_station3Clusters += clusters;
+ break;
+ default:
+ ATH_MSG_FATAL("Unknown tracker station number " << station);
+ break;
+ }
+ }
+
+ // loop over spacepoints and store each space point position
+ SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
+ ATH_CHECK(spacePointContainer.isValid());
+ for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+ m_nspacepoints += spacePointCollection->size();
+ for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) {
+ auto pos = spacePoint->globalPosition();
+ m_spacepointX.push_back(pos.x());
+ m_spacepointY.push_back(pos.y());
+ m_spacepointZ.push_back(pos.z());
+ }
+ }
+
+ // loop over track segments and store position, momentum, chi2, and nDOF for each segment
+ SG::ReadHandle<TrackCollection> trackSegmentCollection {m_trackSegmentCollection, ctx};
+ ATH_CHECK(trackSegmentCollection.isValid());
+ for (const Trk::Track* trackSeg : *trackSegmentCollection) {
+ if (trackSeg == nullptr) continue;
+ m_ntracksegs += 1;
+ m_trackseg_Chi2.push_back(trackSeg->fitQuality()->chiSquared());
+ m_trackseg_DoF.push_back(trackSeg->fitQuality()->numberDoF());
+ auto SegParameters = trackSeg->trackParameters()->front();
+ const Amg::Vector3D SegPosition = SegParameters->position();
+ const Amg::Vector3D SegMomentum = SegParameters->momentum();
+ m_trackseg_x.push_back(SegPosition.x());
+ m_trackseg_y.push_back(SegPosition.y());
+ m_trackseg_z.push_back(SegPosition.z());
+ m_trackseg_px.push_back(SegMomentum.x());
+ m_trackseg_py.push_back(SegMomentum.y());
+ m_trackseg_pz.push_back(SegMomentum.z());
+ }
+
// finished processing event, now fill ntuple tree
m_tree->Fill();
m_eventsPassed += 1;
return StatusCode::SUCCESS;
}
-
StatusCode NtupleDumperAlg::finalize()
{
- ATH_MSG_INFO("Number of events passed Ntuple selectioon = " << m_eventsPassed);
+ ATH_MSG_INFO("Number of events passed Ntuple selection = " << m_eventsPassed);
+ ATH_MSG_INFO("Number of events failing GRL selection = " << m_eventsFailedGRL);
return StatusCode::SUCCESS;
}
-bool NtupleDumperAlg::waveformHitOK(const xAOD::WaveformHit* hit) const
-{
- if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED) || hit->status_bit(xAOD::WaveformStatus::SECONDARY)) return false;
- return true;
-}
-
void
NtupleDumperAlg::clearTree() const
{
- // set all float variables to NaN
- // set all int variables to 999999 (can't set int to NaN, because NaN is a double)
- // set all counter variables to zero
- // set all trigger words to zero
- m_run_number = 999999;
- m_event_number = 999999;
- m_event_time = 999999;
- m_bcid = 999999;
-
- m_fillNumber = 999999;
- m_betaStar = NaN;
- m_crossingAngle = NaN;
+ // don't use NaN, they are annoying when doing cuts
+ m_run_number = 0;
+ m_event_number = 0;
+ m_event_time = 0;
+ m_bcid = 0;
+ m_in_grl = 0;
+
+ m_fillNumber = 0;
+ m_betaStar = 0;
+ m_crossingAngle = 0;
m_distanceToCollidingBCID = 999999;
m_distanceToUnpairedB1 = 999999;
m_distanceToUnpairedB2 = 999999;
@@ -1156,25 +1266,30 @@ NtupleDumperAlg::clearTree() const
m_inputBitsNext=0;
for(int ii=0;ii<15;ii++) {
- m_wave_localtime[ii]=NaN;
- m_wave_peak[ii]=NaN;
- m_wave_width[ii]=NaN;
- m_wave_charge[ii]=NaN;
-
- m_wave_raw_peak[ii]=NaN;
- m_wave_raw_charge[ii]=NaN;
- m_wave_baseline_mean[ii]=NaN;
- m_wave_baseline_rms[ii]=NaN;
+ m_wave_localtime[ii]=0;
+ m_wave_peak[ii]=0;
+ m_wave_width[ii]=0;
+ m_wave_charge[ii]=0;
+
+ m_wave_raw_peak[ii]=0;
+ m_wave_raw_charge[ii]=0;
+ m_wave_baseline_mean[ii]=0;
+ m_wave_baseline_rms[ii]=0;
m_wave_status[ii]=1; // default = 1 means below threshold
- m_calibrated_nMIP[ii]=NaN;
- m_calibrated_E_dep[ii]=NaN;
- m_calibrated_E_EM[ii]=NaN;
+ m_calibrated_nMIP[ii]=0;
+ m_calibrated_E_dep[ii]=0;
+ m_calibrated_E_EM[ii]=0;
}
+ m_scintHit = 0;
+
m_calo_total_nMIP=0;
m_calo_total_E_dep=0;
m_calo_total_E_EM=0;
+ m_calo_total_raw_E_EM=0;
+ m_calo_total_E_EM_fudged=0;
+ m_calo_total_raw_E_EM_fudged=0;
m_preshower_total_nMIP=0;
m_preshower_total_E_dep=0;
@@ -1264,6 +1379,11 @@ NtupleDumperAlg::clearTree() const
m_thetaxCalo.clear();
m_thetayCalo.clear();
+ m_x_atMaxRadius.clear();
+ m_y_atMaxRadius.clear();
+ m_z_atMaxRadius.clear();
+ m_r_atMaxRadius.clear();
+
m_t_pdg.clear();
m_t_barcode.clear();
m_t_truthHitRatio.clear();
@@ -1292,7 +1412,7 @@ NtupleDumperAlg::clearTree() const
m_truthParticleIsFiducial.clear();
m_truthLeptonMomentum = 0;
- m_truthBarcode = 0;
+ m_truthBarcode = -1;
m_truthPdg = 0;
m_truth_P.clear();
@@ -1334,28 +1454,32 @@ NtupleDumperAlg::clearTree() const
}
-void NtupleDumperAlg::setNaN() const {
+void NtupleDumperAlg::clearTrackTruth() const {
m_t_pdg.push_back(0);
m_t_barcode.push_back(-1);
- m_t_truthHitRatio.push_back(NaN);
- m_t_prodVtx_x.push_back(NaN);
- m_t_prodVtx_y.push_back(NaN);
- m_t_prodVtx_z.push_back(NaN);
- m_t_decayVtx_x.push_back(NaN);
- m_t_decayVtx_y.push_back(NaN);
- m_t_decayVtx_z.push_back(NaN);
- m_t_px.push_back(NaN);
- m_t_py.push_back(NaN);
- m_t_pz.push_back(NaN);
- m_t_theta.push_back(NaN);
- m_t_phi.push_back(NaN);
- m_t_p.push_back(NaN);
- m_t_pT.push_back(NaN);
- m_t_eta.push_back(NaN);
+ m_t_truthHitRatio.push_back(0);
+ m_t_prodVtx_x.push_back(999999);
+ m_t_prodVtx_y.push_back(999999);
+ m_t_prodVtx_z.push_back(999999);
+ m_t_decayVtx_x.push_back(999999);
+ m_t_decayVtx_y.push_back(999999);
+ m_t_decayVtx_z.push_back(999999);
+ m_t_px.push_back(0);
+ m_t_py.push_back(0);
+ m_t_pz.push_back(0);
+ m_t_theta.push_back(999999);
+ m_t_phi.push_back(999999);
+ m_t_p.push_back(0);
+ m_t_pT.push_back(0);
+ m_t_eta.push_back(999999);
for (int station = 0; station < 4; ++station) {
- m_t_st_x[station].push_back(NaN);
- m_t_st_y[station].push_back(NaN);
- m_t_st_z[station].push_back(NaN);
+ m_t_st_x[station].push_back(999999);
+ m_t_st_y[station].push_back(999999);
+ m_t_st_z[station].push_back(999999);
}
m_isFiducial.push_back(false);
}
+
+double NtupleDumperAlg::radius(const Acts::Vector3 &position) const {
+ return std::sqrt(position.x() * position.x() + position.y() * position.y());
+}
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 6d4ccc381a4b65831678a657c035ecc4e79c2b65..339c0838b22047dd146d71163acebdc472f7d1f3 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -25,6 +25,7 @@
#include "StoreGate/ReadDecorHandle.h"
#include <vector>
+#include <nlohmann/json.hpp>
class TTree;
class TH1;
@@ -52,12 +53,13 @@ private:
bool waveformHitOK(const xAOD::WaveformHit* hit) const;
void clearTree() const;
- void setNaN() const;
+ void clearTrackTruth() const;
void addBranch(const std::string &name,float* var);
void addBranch(const std::string &name,unsigned int* var);
void addWaveBranches(const std::string &name, int nchannels, int first);
void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
void addCalibratedBranches(const std::string &name, int nchannels, int first);
+ double radius(const Acts::Vector3 &position) const;
ServiceHandle <ITHistSvc> m_histSvc;
@@ -101,17 +103,29 @@ 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" };
- IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
- DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
- 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"};
+ 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 no veto signal adn Calo signal > 20 GeV e-" };
+ BooleanProperty m_onlyBlinded { this, "OnlyBlinded", false, "Only events that would be blinded are saved" };
+ BooleanProperty m_doTrigFilter { this, "DoTrigFilter", false, "Only events that pass trigger cuts are passed" };
+ BooleanProperty m_doScintFilter { this, "DoScintFilter", false, "Only events that pass scintillator coincidence cuts are passed" };
+ BooleanProperty m_doTrackFilter { this, "DoTrackFilter", true, "Only events that have >= 1 long track are passed, also non-colliding events with a track or calo signal are passed" };
- DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"};
+ 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" };
+ IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
+ DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
+ 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", 100000.0, "Blind events with Ecal energy above threshold (in MeV)"};
+ DoubleProperty m_caloMC_FudgeFactor {this, "CaloFudgeFactorMC", 1.088, "Correct the MC energy calibration by this fudge factor"};
+
+ BooleanProperty m_applyGoodRunsList {this, "ApplyGoodRunsList", false, "Only write out events passing GRL (data only)"};
+ StringProperty m_goodRunsList {this, "GoodRunsList", "", "GoodRunsList in json format"};
+
+ // json object to hold data read from GRL file (or empty if not)
+ nlohmann::json m_grl;
double m_baseEventCrossSection {1.0};
const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -124,6 +138,7 @@ private:
mutable unsigned int m_event_number;
mutable unsigned int m_event_time;
mutable unsigned int m_bcid;
+ mutable unsigned int m_in_grl;
mutable unsigned int m_fillNumber;
mutable float m_betaStar;
@@ -150,7 +165,9 @@ private:
mutable float m_wave_baseline_mean[15];
mutable float m_wave_baseline_rms[15];
mutable unsigned int m_wave_status[15];
-
+
+ mutable unsigned int m_scintHit;
+
mutable float m_calibrated_nMIP[15];
mutable float m_calibrated_E_dep[15];
mutable float m_calibrated_E_EM[15];
@@ -158,6 +175,10 @@ private:
mutable float m_calo_total_nMIP;
mutable float m_calo_total_E_dep;
mutable float m_calo_total_E_EM;
+ mutable float m_calo_total_raw_E_EM;
+
+ mutable float m_calo_total_E_EM_fudged;
+ mutable float m_calo_total_raw_E_EM_fudged;
mutable float m_preshower_total_nMIP;
mutable float m_preshower_total_E_dep;
@@ -170,9 +191,6 @@ private:
mutable float m_Preshower12_Edep;
mutable float m_Preshower13_Edep;
- mutable float m_MIP_sim_Edep_calo;
- mutable float m_MIP_sim_Edep_preshower;
-
mutable float m_clock_phase;
mutable unsigned int m_station0Clusters;
@@ -247,6 +265,10 @@ private:
mutable std::vector<double> m_yCalo;
mutable std::vector<double> m_thetaxCalo;
mutable std::vector<double> m_thetayCalo;
+ mutable std::vector<double> m_x_atMaxRadius;
+ mutable std::vector<double> m_y_atMaxRadius;
+ mutable std::vector<double> m_z_atMaxRadius;
+ mutable std::vector<double> m_r_atMaxRadius;
mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
@@ -324,16 +346,13 @@ private:
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;
mutable double m_crossSection;
mutable int m_eventsPassed = 0;
-
+ mutable int m_eventsFailedGRL = 0;
};
inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {
diff --git a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
index d7ea3c3ab921b57e3df3c5c4be56c102f873c8f7..5ff7e6b96d1c04ce7869062feb8bea942e8f380e 100755
--- a/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/G4FaserAlgConfigNew_Test.py
@@ -20,8 +20,8 @@ if __name__ == '__main__':
# Import and set config flags
#
from CalypsoConfiguration.AllConfigFlags import ConfigFlags
- ConfigFlags.Exec.MaxEvents = 4 # can be overridden from command line with --maxEvt=<number>
- ConfigFlags.Exec.SkipEvents = 0 # can be overridden from command line with --skipEvt=<number>
+ ConfigFlags.Exec.MaxEvents = 4 # can be overridden from command line with --evtMax=<number>
+ ConfigFlags.Exec.SkipEvents = 0 # can be overridden from command line with --skipEvents=<number>
from AthenaConfiguration.Enums import ProductionStep
ConfigFlags.Common.ProductionStep = ProductionStep.Simulation
#
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8d8665ba924345e2fdcdb409384cdb66291446c3
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(OverlayRDO)
+
+atlas_add_component(
+ OverlayRDO
+ src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h src/HistoRDOAlg.h
+ src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx src/HistoRDOAlg.cxx
+ src/component/OverlayRDO_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfb2dfefe48927be61f6416bd206a78e0b754eab
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
@@ -0,0 +1,127 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def HistoRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ HistoRDOAlg = CompFactory.HistoRDOAlg("HistoRDOAlg",**kwargs)
+ acc.addEventAlgo(HistoRDOAlg)
+ HistoRDOAlg.OutputLevel = VERBOSE
+ # ItemList = []
+ # # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ # ItemList += ["xAOD::EventInfo#EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ # ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+ # ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+ # ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+ # # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ # acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST2 DATAFILE='RDOtree.root' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ # ConfigFlags.Overlay.DataOverlay = False
+ ConfigFlags.Input.Files = [
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin0RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin1RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin2RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin3RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin4RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin5RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin6RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin7RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin8RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin0RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin1RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin2RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin3RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin4RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin5RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin6RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin7RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin8RDO.pool.root',
+ ]
+ # ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+ # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+ # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+ # ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(HistoRDOAlgCfg(ConfigFlags))
+
+ # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ # acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+ # ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..bce33065b965395662707d54a82606d1e7038b9e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -0,0 +1,106 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def OverlayRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ OverlayRDOAlg = CompFactory.OverlayRDOAlg("OverlayRDOAlg",**kwargs)
+ acc.addEventAlgo(OverlayRDOAlg)
+
+ ItemList = []
+ # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ ItemList += ["xAOD::EventInfo#EventInfo"]
+ ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+ ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+ ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+ # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+
+
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Overlay.DataOverlay = False
+ ConfigFlags.Input.Files = [ 'Pos_RDO.pool.root']
+ ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+ # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+ # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+ ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(OverlayRDOAlgCfg(ConfigFlags))
+
+ # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ # acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+ # ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..29534a9b3293511b430305770da134c073b10213
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -0,0 +1,156 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+chargePrefix = "Pos_"
+
+def SelectRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ SelectRDOAlg = CompFactory.SelectRDOAlg("SelectRDOAlg",**kwargs)
+ SelectRDOAlg.AcceptPositive = ("Pos_" in chargePrefix)
+ acc.addEventAlgo(SelectRDOAlg)
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin0RDO",Xmin=-66, Xmax=0, Ymin=-66, Ymax=0))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin1RDO",Xmin=0, Xmax=66, Ymin=-66, Ymax=0))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin2RDO",Xmin=-66, Xmax=0, Ymin=0, Ymax=66))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin3RDO",Xmin=0, Xmax=66, Ymin=0, Ymax=66))
+
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin4RDO",Xmin=-33, Xmax=33, Ymin=-33, Ymax=33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin5RDO",Xmin=-33, Xmax=33, Ymin=33, Ymax=99))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin6RDO",Xmin=-33, Xmax=33, Ymin=-99, Ymax=-33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin7RDO",Xmin=-99, Xmax=-33, Ymin=-33, Ymax=33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin8RDO",Xmin=33, Xmax=99, Ymin=-33, Ymax=33))
+
+ ItemList = []
+ # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ ItemList += ["xAOD::EventInfo#EventInfo"]
+ ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+ ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_RDOs"]
+ ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_EDGEMODE_RDOs"]
+ # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin0RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo0 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin0RDO")
+ osrdo0.RequireAlgs += ["SelectRDOAlg","Bin0RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin1RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo1 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin1RDO")
+ osrdo1.RequireAlgs += ["SelectRDOAlg","Bin1RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin2RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo2 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin2RDO")
+ osrdo2.RequireAlgs += ["SelectRDOAlg","Bin2RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin3RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo3 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin3RDO")
+ osrdo3.RequireAlgs += ["SelectRDOAlg","Bin3RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin4RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo4 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin4RDO")
+ osrdo4.RequireAlgs += ["SelectRDOAlg","Bin4RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin5RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo5 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin5RDO")
+ osrdo5.RequireAlgs += ["SelectRDOAlg","Bin5RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin6RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo6 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin6RDO")
+ osrdo6.RequireAlgs += ["SelectRDOAlg","Bin6RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin7RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo7 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin7RDO")
+ osrdo7.RequireAlgs += ["SelectRDOAlg","Bin7RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin8RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo8 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin8RDO")
+ osrdo8.RequireAlgs += ["SelectRDOAlg","Bin8RDO"]
+
+ return acc
+
+
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [ #'/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root'
+ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root'
+ ]
+ # ConfigFlags.Input.Files = [ 'Faser-Physics-009171-00006-r0013-xAOD.root',
+ # 'Faser-Physics-009166-00485-r0013-xAOD.root']
+ # ConfigFlags.Output.RDOFileName = chargePrefix + "RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(SelectRDOAlgCfg(ConfigFlags, TrackCollection = "CKFTrackCollectionWithoutIFT", OutputTrackCollection = chargePrefix + "CKFTrackCollectionWithoutIFT"))
+
+ from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + chargePrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + chargePrefix + "EventInfoAux.",
+ # "TrackCollection#CKFTrackCollectionWithoutIFT->" + chargePrefix + "CKFTrackCollectionWithoutIFT",
+ "FaserSCT_RDO_Container#SCT_RDOs->" + chargePrefix + "SCT_RDOs",
+ "FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs->" + chargePrefix + "SCT_EDGEMODE_RDOs",
+ # "Tracker::FaserSCT_ClusterContainer#SCT_ClusterContainer->" + chargePrefix + "SCT_ClusterContainer"
+ ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a917447dde480b9b82a316653fb8e76cf76a2156
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
@@ -0,0 +1,55 @@
+#include "BinRDOAlg.h"
+#include "TrkTrack/Track.h"
+
+BinRDOAlg::BinRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode BinRDOAlg::initialize()
+{
+ ATH_CHECK(m_trackCollection.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode BinRDOAlg::execute(const EventContext &ctx) const
+{
+
+ setFilterPassed(false, ctx);
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+ // track cuts
+
+ double xUpstream {0.0};
+ double yUpstream {0.0};
+
+ for (const Trk::Track* track : *trackCollection)
+ {
+ if (track == nullptr) continue;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ }
+ xUpstream = upstreamParameters->position().x();
+ yUpstream = upstreamParameters->position().y();
+ break;
+ }
+
+
+ // set filter passed
+
+ setFilterPassed(((xUpstream >= m_xMin) && (xUpstream <= m_xMax) && (yUpstream >= m_yMin) && (yUpstream <= m_yMax)), ctx);
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode BinRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8e86e844adf11e7a1a4021b9e8d7a11757ce58c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
@@ -0,0 +1,23 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+
+class BinRDOAlg : public AthReentrantAlgorithm {
+public:
+ BinRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~BinRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ DoubleProperty m_xMin {this, "Xmin", -25.0, "Minimum x position to accept"};
+ DoubleProperty m_xMax {this, "Xmax", 25.0, "Maximum x position to accept"};
+ DoubleProperty m_yMin {this, "Ymin", -25.0, "Minimum y position to accept"};
+ DoubleProperty m_yMax {this, "Ymax", 25.0, "Maximum y position to accept"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..78096c60004a9e183bf83aa235ef72501ddee4c4
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -0,0 +1,247 @@
+#include "HistoRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include <TTree.h>
+
+
+HistoRDOAlg::HistoRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator),
+ AthHistogramming(name),
+ m_histSvc("THistSvc/THistSvc", name) {}
+
+StatusCode HistoRDOAlg::initialize()
+{
+ m_tree = new TTree("overlay", "Overlay tree");
+
+ m_tree->Branch("longTracks1", &m_longTracks1, "longTracks1/I");
+ m_tree->Branch("x1", &m_x1);
+ m_tree->Branch("y1", &m_y1);
+ m_tree->Branch("p1", &m_p1);
+ m_tree->Branch("theta1", &m_theta1);
+ m_tree->Branch("phi1", &m_phi1);
+ m_tree->Branch("charge1", &m_charge1);
+ m_tree->Branch("chiSquared1", &m_chiSquared1);
+ m_tree->Branch("nDoF1", &m_nDoF1);
+ m_tree->Branch("pPolar1", &m_pTotPolar1, "pPolar1/D");
+
+ m_tree->Branch("longTracks2", &m_longTracks2, "longTracks2/I");
+ m_tree->Branch("x2", &m_x2);
+ m_tree->Branch("y2", &m_y2);
+ m_tree->Branch("p2", &m_p2);
+ m_tree->Branch("theta2", &m_theta2);
+ m_tree->Branch("phi2", &m_phi2);
+ m_tree->Branch("charge2", &m_charge2);
+ m_tree->Branch("chiSquared2", &m_chiSquared2);
+ m_tree->Branch("nDoF2", &m_nDoF2);
+ m_tree->Branch("pPolar2", &m_pTotPolar2, "pPolar2/D");
+
+ m_tree->Branch("yAgreement", &m_yAgreement, "yAgreement/D");
+ m_tree->Branch("matchFraction", &m_matchFraction, "matchFraction/D");
+
+ ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+ ATH_CHECK(m_trackCollection1.initialize());
+ ATH_CHECK(m_trackCollection2.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
+{
+
+ m_longTracks1 = 0;
+ m_x1.clear();
+ m_y1.clear();
+ m_p1.clear();
+ m_theta1.clear();
+ m_phi1.clear();
+ m_charge1.clear();
+ m_chiSquared1.clear();
+ m_nDoF1.clear();
+ m_pTotPolar1 = 0;
+
+ m_longTracks2 = 0;
+ m_x2.clear();
+ m_y2.clear();
+ m_p2.clear();
+ m_theta2.clear();
+ m_phi2.clear();
+ m_charge2.clear();
+ m_chiSquared2.clear();
+ m_nDoF2.clear();
+ m_pTotPolar2 = 0;
+
+ m_yAgreement = 0;
+ m_matchFraction = 0;
+
+ SG::ReadHandle<TrackCollection> trackCollection1 {m_trackCollection1, ctx};
+ // ATH_CHECK(trackCollection1.isValid());
+ SG::ReadHandle<TrackCollection> trackCollection2 {m_trackCollection2, ctx};
+ // ATH_CHECK(trackCollection2.isValid());
+
+ std::vector<const Trk::TrackParameters*> upstreamSingle {};
+ if (trackCollection1.isValid())
+ {
+ Amg::Vector3D pTot {0, 0, 0};
+ for (const Trk::Track* track : *trackCollection1)
+ {
+ if (track == nullptr) continue;
+ m_longTracks1++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ ATH_MSG_VERBOSE("TrackCollection 1, Track " << m_longTracks1 << ", z = " << params->position().z() << ", y = " << params->position().y() << ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+ }
+ upstreamSingle.push_back(upstreamParameters);
+ m_x1.push_back(upstreamParameters->position().x());
+ m_y1.push_back(upstreamParameters->position().y());
+ m_p1.push_back(upstreamParameters->momentum().mag());
+ m_charge1.push_back(upstreamParameters->charge());
+ m_theta1.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+ m_phi1.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+ m_chiSquared1.push_back(track->fitQuality()->chiSquared());
+ m_nDoF1.push_back(track->fitQuality()->doubleNumberDoF());
+ pTot += upstreamParameters->momentum();
+ }
+ m_pTotPolar1 = asin(pTot.perp()/pTot.mag());
+ }
+
+ std::vector<const Trk::TrackParameters*> upstreamOverlay {};
+ if (trackCollection2.isValid())
+ {
+ Amg::Vector3D pTot {0, 0, 0};
+ for (const Trk::Track* track : *trackCollection2)
+ {
+ if (track == nullptr) continue;
+ m_longTracks2++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ ATH_MSG_VERBOSE("TrackCollection 2, Track " << m_longTracks2 << ", z = " << params->position().z()<< ", y = " << params->position().y()<< ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+ }
+ upstreamOverlay.push_back(upstreamParameters);
+ m_x2.push_back(upstreamParameters->position().x());
+ m_y2.push_back(upstreamParameters->position().y());
+ m_p2.push_back(upstreamParameters->momentum().mag());
+ m_charge2.push_back(upstreamParameters->charge());
+ m_theta2.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+ m_phi2.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+ m_chiSquared2.push_back(track->fitQuality()->chiSquared());
+ m_nDoF2.push_back(track->fitQuality()->doubleNumberDoF());
+ pTot += upstreamParameters->momentum();
+ }
+ m_pTotPolar2 = asin(pTot.perp()/pTot.mag());
+ }
+
+ if (trackCollection1.isValid() && trackCollection2.isValid() && m_longTracks2 ==2)
+ {
+ double qOverPSingle1 = upstreamSingle[0]->charge()/upstreamSingle[0]->momentum().mag();
+ double qOverPCovSingle1 = (*upstreamSingle[0]->covariance())(4,4);
+ double qOverPSingle2 = upstreamSingle[1]->charge()/upstreamSingle[1]->momentum().mag();
+ double qOverPCovSingle2 = (*upstreamSingle[1]->covariance())(4,4);
+ double qOverPOverlay1 = upstreamOverlay[0]->charge()/upstreamOverlay[0]->momentum().mag();
+ double qOverPCovOverlay1 = (*upstreamOverlay[0]->covariance())(4,4);
+ double qOverPOverlay2 = upstreamOverlay[1]->charge()/upstreamOverlay[1]->momentum().mag();
+ double qOverPCovOverlay2 = (*upstreamOverlay[1]->covariance())(4,4);
+
+ double agreementDirect = pow(qOverPSingle1 - qOverPOverlay1, 2) / (qOverPCovSingle1 + qOverPCovOverlay1) +
+ pow(qOverPSingle2 - qOverPOverlay2, 2) / (qOverPCovSingle2 + qOverPCovOverlay2);
+ double agreementSwap = pow(qOverPSingle1 - qOverPOverlay2, 2) / (qOverPCovSingle1 + qOverPCovOverlay2) +
+ pow(qOverPSingle2 - qOverPOverlay1, 2) / (qOverPCovSingle2 + qOverPCovOverlay1);
+
+ const Trk::Track* firstSingle {nullptr};
+ const Trk::Track* firstOverlay {nullptr};
+ const Trk::Track* secondSingle {nullptr};
+ const Trk::Track* secondOverlay {nullptr};
+ if (agreementDirect <= agreementSwap)
+ {
+ firstSingle = (*trackCollection1)[0];
+ firstOverlay = (*trackCollection2)[0];
+ secondSingle = (*trackCollection1)[1];
+ secondOverlay = (*trackCollection2)[1];
+ ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay1 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay2);
+ }
+ else
+ {
+ firstSingle = (*trackCollection1)[0];
+ firstOverlay = (*trackCollection2)[1];
+ secondSingle = (*trackCollection1)[1];
+ secondOverlay = (*trackCollection2)[0];
+ ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay2 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay1);
+ }
+ int nZMatch {0};
+ auto singleParameters = firstSingle->trackParameters();
+ auto overlayParameters = firstOverlay->trackParameters();
+ for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+ {
+ auto singleState = (*singleParameters)[iSingle];
+ auto overlayState = (*overlayParameters)[iOverlay];
+ auto singleCov = singleState->covariance();
+ auto overlayCov = overlayState->covariance();
+ if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+ {
+ nZMatch++;
+ m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+ iSingle++;
+ iOverlay++;
+ ATH_MSG_VERBOSE("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
+ }
+ else if (singleState->position().z() < overlayState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iSingle++;
+ }
+ else if (overlayState->position().z() < singleState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iOverlay++;
+ }
+ }
+ double matchFract1 { ((double) nZMatch)/singleParameters->size() };
+ ATH_MSG_VERBOSE("MatchFract1: " << matchFract1);
+ nZMatch = 0;
+ singleParameters = secondSingle->trackParameters();
+ overlayParameters = secondOverlay->trackParameters();
+ for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+ {
+ auto singleState = (*singleParameters)[iSingle];
+ auto overlayState = (*overlayParameters)[iOverlay];
+ auto singleCov = singleState->covariance();
+ auto overlayCov = overlayState->covariance();
+ if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+ {
+ nZMatch++;
+ m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+ iSingle++;
+ iOverlay++;
+ ATH_MSG_VERBOSE("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
+ }
+ else if (singleState->position().z() < overlayState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iSingle++;
+ }
+ else if (overlayState->position().z() < singleState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iOverlay++;
+ }
+ }
+ double matchFract2 { ((double) nZMatch)/singleParameters->size() };
+ ATH_MSG_VERBOSE("MatchFract2: " << matchFract2);
+ m_matchFraction = std::min(matchFract1, matchFract2);
+ ATH_MSG_VERBOSE("MatchFraction: " << m_matchFraction);
+ }
+
+ m_tree->Fill();
+ return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..4fe18f7ad71e075b8b37b56518058da6b7be6f0e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
@@ -0,0 +1,58 @@
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+
+#include <vector>
+
+class TTree;
+
+class HistoRDOAlg : public AthReentrantAlgorithm, AthHistogramming {
+public:
+ HistoRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~HistoRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+
+ ServiceHandle <ITHistSvc> m_histSvc;
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection1 { this, "TrackCollection1", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection 1 name" };
+ SG::ReadHandleKey<TrackCollection> m_trackCollection2 { this, "TrackCollection2", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection 2 name" };
+
+// mutable TTree* m_tree;
+
+ mutable TTree* m_tree;
+
+ mutable int m_longTracks1;
+ mutable std::vector<double> m_p1;
+ mutable std::vector<double> m_theta1;
+ mutable std::vector<double> m_phi1;
+ mutable std::vector<double> m_x1;
+ mutable std::vector<double> m_y1;
+ mutable std::vector<double> m_charge1;
+ mutable std::vector<double> m_chiSquared1;
+ mutable std::vector<double> m_nDoF1;
+ mutable double m_pTotPolar1;
+
+ mutable int m_longTracks2;
+ mutable std::vector<double> m_p2;
+ mutable std::vector<double> m_theta2;
+ mutable std::vector<double> m_phi2;
+ mutable std::vector<double> m_x2;
+ mutable std::vector<double> m_y2;
+ mutable std::vector<double> m_charge2;
+ mutable std::vector<double> m_chiSquared2;
+ mutable std::vector<double> m_nDoF2;
+ mutable double m_pTotPolar2;
+
+ mutable double m_yAgreement;
+ mutable double m_matchFraction;
+
+};
+
+inline const ServiceHandle <ITHistSvc> &HistoRDOAlg::histSvc() const {
+ return m_histSvc;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9b88a96f9590462d5c03d6ffba601598185ec109
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -0,0 +1,186 @@
+#include "OverlayRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include <cmath>
+
+
+
+OverlayRDOAlg::OverlayRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+ {}
+
+StatusCode OverlayRDOAlg::initialize()
+{
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+ ATH_CHECK(m_posTrackKey.initialize());
+ ATH_CHECK(m_negTrackKey.initialize());
+ ATH_CHECK(m_posRdoKey.initialize());
+ ATH_CHECK(m_negRdoKey.initialize());
+ ATH_CHECK(m_posEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_negEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_outRdoKey.initialize());
+ ATH_CHECK(m_outEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
+{
+
+ SG::ReadHandle<TrackCollection> posTrackCollection {m_posTrackKey, ctx};
+ ATH_CHECK(posTrackCollection.isValid());
+ SG::ReadHandle<TrackCollection> negTrackCollection {m_negTrackKey, ctx};
+ ATH_CHECK(negTrackCollection.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> posRdoContainer {m_posRdoKey, ctx};
+ ATH_CHECK(posRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> negRdoContainer {m_negRdoKey, ctx};
+ ATH_CHECK(negRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> posEdgeModeRdoContainer {m_posEdgeModeRdoKey, ctx};
+ ATH_CHECK(posEdgeModeRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> negEdgeModeRdoContainer {m_negEdgeModeRdoKey, ctx};
+ ATH_CHECK(negEdgeModeRdoContainer.isValid());
+
+ // SCT_RDOs
+
+ std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> > rdoDB {};
+ processContainer(*posRdoContainer, rdoDB);
+ processContainer(*negRdoContainer, rdoDB);
+
+ SG::WriteHandle<FaserSCT_RDO_Container> outRdoContainer {m_outRdoKey, ctx};
+ ATH_CHECK(outRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+ processDB(rdoDB, outRdoContainer, [](int i) {return true;});
+
+ // SCT_EDGEMODE_RDOs
+
+ rdoDB.clear();
+ processContainer(*posEdgeModeRdoContainer, rdoDB);
+ processContainer(*negEdgeModeRdoContainer, rdoDB);
+
+ SG::WriteHandle<FaserSCT_RDO_Container> outEdgeModeRdoContainer {m_outEdgeModeRdoKey, ctx};
+ ATH_CHECK(outEdgeModeRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+ processDB(rdoDB, outEdgeModeRdoContainer, [](int hitPattern) { return (((hitPattern & 0x2) == 0 ) || ((hitPattern & 0x4) != 0) ) ? false : true;});
+
+ // Tracks
+
+ SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ for (auto theTrack : *posTrackCollection)
+ {
+ outputTracks->push_back(cloneTrack(theTrack) );
+ }
+
+ for (auto theTrack : *negTrackCollection)
+ {
+ outputTracks->push_back(cloneTrack(theTrack) );
+ }
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+ return StatusCode::SUCCESS;
+}
+
+
+void
+OverlayRDOAlg::processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB, SG::WriteHandle<FaserSCT_RDO_Container>& outRdoContainer, std::function<bool(int)> lambda) const
+{
+ for (auto& entry : rdoDB)
+ {
+ IdentifierHash waferHash = entry.first;
+ Identifier id = m_sctHelper->wafer_id(waferHash);
+ auto vec = entry.second;
+ if (vec.size() == 0)
+ {
+ std::cout << "Unexpected zero-length collection vector for ID = " << id << " ..." << std::endl;
+ continue;
+ }
+ std::unique_ptr<FaserSCT_RDO_Collection> current_collection = std::make_unique<FaserSCT_RDO_Collection>(waferHash);
+ current_collection->setIdentifier(id);
+ std::map<int, unsigned int> stripMap;
+ for (auto& collection : vec)
+ {
+ for (auto rawData : *collection)
+ {
+ // const FaserSCT3_RawData* data = dynamic_cast<const FaserSCT3_RawData*>(rawData);
+ Identifier stripID = rawData->identify();
+ int stripNumber = m_sctHelper->strip(stripID);
+ int groupSize = rawData->getGroupSize();
+ unsigned int dataWord = rawData->getWord();
+ // int time = data->getTimeBin();
+ // std::cout << "Word: " << std::hex << dataWord << std::dec << " Time: " << time << std::endl;
+ for (int i = stripNumber; i < stripNumber + groupSize; i++)
+ stripMap[i] |= dataWord;
+ }
+ }
+
+ for (auto stripEntry : stripMap)
+ {
+ Identifier rdoID {m_sctHelper->strip_id(id, stripEntry.first)};
+ if (!lambda(stripEntry.second>>22)) continue;
+ current_collection->emplace_back(new FaserSCT3_RawData(rdoID, stripEntry.second, std::vector<int> {}));
+ }
+
+ outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
+ }
+}
+
+
+Trk::Track* OverlayRDOAlg::cloneTrack(const Trk::Track* theTrack) const
+{
+ Trk::TrackInfo i { theTrack->info() };
+ Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+ DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
+
+ for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+ {
+ const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+ std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+ Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+ s->push_back(newState);
+ }
+ std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+ return new Trk::Track {i, std::move(*sink) , q } ;
+}
+
+void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const
+{
+ size_t nLong{0};
+ size_t nData{0};
+ for( const auto& collection : container)
+ {
+ if (collection->empty()) continue;
+ Identifier id = collection->identify();
+ IdentifierHash hash = m_sctHelper->wafer_hash(id);
+ if (rdoDB.count(hash))
+ {
+ rdoDB[hash].push_back(collection);
+ }
+ else
+ {
+ rdoDB[hash] = std::vector<const FaserSCT_RDO_Collection*>{collection};
+ }
+ for(const auto& rawdata : *collection)
+ {
+ nData++;
+ if (rawdata->getGroupSize() > 1) nLong++;
+ }
+ }
+ std::cout << nLong << " long data out of " << nData << " total RDOs" << std::endl;
+}
+
+
+StatusCode OverlayRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..eba2f90d7cde46873bbe9124e17876fef5059e5a
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -0,0 +1,38 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+
+class TTree;
+class FaserSCT_ID;
+
+class OverlayRDOAlg : public AthReentrantAlgorithm {
+public:
+ OverlayRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~OverlayRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ void processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const;
+ void processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& db, SG::WriteHandle<FaserSCT_RDO_Container>& outContainer, std::function<bool(int)> lambda) const;
+ Trk::Track* cloneTrack(const Trk::Track* originalTrack) const;
+
+ SG::ReadHandleKey<TrackCollection> m_posTrackKey { this, "PosTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_negTrackKey { this, "NegTrackCollection", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_posRdoKey { this, "PosRdoContainer", "Pos_SCT_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_posEdgeModeRdoKey { this, "PosEdgeModeRdoContainer", "Pos_SCT_EDGEMODE_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_negEdgeModeRdoKey { this, "NegEdgeModeRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
+
+ SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputRDOObjectName", "SCT_RDOs", "Output RDO Object name"};
+ SG::WriteHandleKey<FaserSCT_RDO_Container> m_outEdgeModeRdoKey{this, "OutputEDGEMODEObjectName", "SCT_EDGEMODE_RDOs", "Output EDGEMODE Object name"};
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Orig_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+ const FaserSCT_ID* m_sctHelper;
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b33bed9d18ac4d51f2e2dc7ae9fc52217cbffc67
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -0,0 +1,187 @@
+#include "SelectRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include <cmath>
+
+SelectRDOAlg::SelectRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode SelectRDOAlg::initialize()
+{
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_segmentCollection.initialize());
+ ATH_CHECK(m_clusterContainer.initialize());
+ ATH_CHECK(m_triggerContainer.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
+{
+
+ setFilterPassed(false, ctx);
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<TrackCollection> segmentCollection {m_segmentCollection, ctx};
+ if (!segmentCollection.isValid()) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+ if (!clusterContainer.isValid()) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ if (!triggerContainer.isValid()) return StatusCode::SUCCESS;
+
+ // segment cuts
+
+ std::vector<int> segmentCount {0, 0, 0, 0};
+ int nSegments {0};
+
+ for (const Trk::Track* segment : *segmentCollection)
+ {
+ if (segment == nullptr) continue;
+ nSegments++;
+ for (auto measurement : *(segment->measurementsOnTrack()))
+ {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr)
+ {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ if (station >= 0 && station <= 3)
+ {
+ segmentCount[station]++;
+ break;
+ }
+ }
+ }
+ }
+ if ((nSegments > m_maxSegmentsTotal) || (segmentCount[1] > m_maxSegmentsStation) || (segmentCount[2] > m_maxSegmentsStation) || (segmentCount[3] > m_maxSegmentsStation))
+ return StatusCode::SUCCESS;
+
+ // track cuts
+
+ int nTracks {0};
+ int nLongTracks {0};
+ double chi2PerDoF {0};
+ double charge {0};
+ double maxRadius {0};
+ double momentum {0};
+ const Trk::Track* theTrack {nullptr};
+
+ for (const Trk::Track* track : *trackCollection)
+ {
+ if (track == nullptr) continue;
+ nTracks++;
+ std::set<int> stationMap;
+ std::set<std::pair<int, int>> layerMap;
+
+ // Check for hit in the three downstream stations
+ for (auto measurement : *(track->measurementsOnTrack()))
+ {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr)
+ {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ int layer = m_sctHelper->layer(id);
+ stationMap.emplace(station);
+ layerMap.emplace(station, layer);
+ }
+ }
+ if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+ int nLayers = std::count_if(layerMap.begin(), layerMap.end(), [](std::pair<int,int> p){return p.first != 0;});
+ if (nLayers < m_minLayers) continue;
+ nLongTracks++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ double radius = sqrt(pow(params->position().x(), 2) + pow(params->position().y(), 2));
+ if (radius > maxRadius) maxRadius = radius;
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ }
+
+ momentum = upstreamParameters->momentum().mag()/1000;
+ chi2PerDoF = track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF();
+ charge = upstreamParameters->charge();
+ theTrack = track;
+ }
+
+ if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (maxRadius > m_maxRadius))
+ return StatusCode::SUCCESS;
+
+ // cluster cuts
+
+ // std::vector<int> clusterCount {0, 0, 0, 0};
+ // for (auto collection : *clusterContainer)
+ // {
+ // Identifier id = collection->identify();
+ // int station = m_sctHelper->station(id);
+ // if (station >= 0 && station <= 3) clusterCount[station] += collection->size();
+ // }
+ // if ((clusterCount[1] > m_maxClustersStation) || (clusterCount[2] > m_maxClustersStation) || (clusterCount[3] > m_maxClustersStation))
+ // return StatusCode::SUCCESS;
+
+ // waveform charge cuts
+
+ std::vector<double> timingCharge {0, 0, 0, 0};
+ const int baseChannel = 8;
+
+ for (auto hit : *triggerContainer) {
+ if ((hit->hit_status()&2)==0) { // dont store secondary hits as they can overwrite the primary hit
+ int ch=hit->channel();
+ timingCharge[ch - baseChannel] += hit->integral()/50;
+ }
+ }
+ if ((timingCharge[0] + timingCharge[1] > m_maxTimingCharge) || (timingCharge[2] + timingCharge[3] > m_maxTimingCharge))
+ return StatusCode::SUCCESS;
+
+ // passes all cuts
+
+ // copy track without associated clusters
+
+ SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ Trk::TrackInfo i { theTrack->info() };
+ Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+ DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
+
+ for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+ {
+ const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+ std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+ Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+ s->push_back(newState);
+ }
+ std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+ outputTracks->push_back(new Trk::Track {i, std::move(*sink) , q } );
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+ // set filter passed
+
+ setFilterPassed(true, ctx);
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..842282d534a23159625a1dabf59c3fbc0494c343
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -0,0 +1,43 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHit.h"
+#include <xAODEventInfo/EventInfo.h>
+
+
+class FaserSCT_ID;
+
+class SelectRDOAlg : public AthReentrantAlgorithm {
+public:
+ SelectRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~SelectRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ const FaserSCT_ID* m_sctHelper;
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_segmentCollection { this, "SegmentCollection", "SegmentFit", "Track segment collection name" };
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+ IntegerProperty m_minLayers {this, "MinLayers", 7, "Minimum hit layers to accept track"};
+ DoubleProperty m_maxChi2PerDoF {this, "MaxChi2PerDoF", 25.0, "Maximum chi2 per degree of freedom to accept track"};
+ DoubleProperty m_minMomentum {this, "MinMomentumGeV", 50.0, "Minimum momentum in GeV to accept track"};
+ DoubleProperty m_maxMomentum {this, "MaxMomentumGeV", 5000.0, "Maximum momentum in GeV to accept track"};
+ DoubleProperty m_maxRadius {this, "MaxRadiusMm", 95.0, "Maximum radius at first measurement to accept track"};
+ IntegerProperty m_maxSegmentsTotal {this, "MaxSegmentsTotal", 4, "Maximum number of segments in three stations to accept track"};
+ IntegerProperty m_maxSegmentsStation {this, "MaxSegmentsStation", 2, "Maximum number of segments in any single station to accept track"};
+ // IntegerProperty m_maxClustersStation {this, "MaxClustersStation", 9, "Maximum number of clusters in any single station to accept track"};
+ DoubleProperty m_maxTimingCharge {this, "MaxTimingChargePc", 70.0, "Maximum charge in pC recorded by upper or lower timing scintillator"};
+ BooleanProperty m_acceptPositive {this, "AcceptPositive", true, "Accept positive (true) or negative (false) charged tracks"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..13e61cb84654d0c3e4096ffa01b0f6efbeb6184c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -0,0 +1,9 @@
+#include "../OverlayRDOAlg.h"
+#include "../SelectRDOAlg.h"
+#include "../BinRDOAlg.h"
+#include "../HistoRDOAlg.h"
+
+DECLARE_COMPONENT(OverlayRDOAlg)
+DECLARE_COMPONENT(SelectRDOAlg)
+DECLARE_COMPONENT(BinRDOAlg)
+DECLARE_COMPONENT(HistoRDOAlg)
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
index f0de1582f2ccf47270203d4d1d6eb9a6bca341e9..eeaadd061bbc5b1c3c14af0a1e128e529d687cae 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
@@ -15,9 +15,9 @@ def SegmentFitAlgBasicCfg(flags, **kwargs):
Tracker__SegmentFitAlg=CompFactory.Tracker.SegmentFitAlg
acc.addEventAlgo(Tracker__SegmentFitAlg(**kwargs))
- thistSvc = CompFactory.THistSvc()
- thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
- acc.addService(thistSvc)
+ # thistSvc = CompFactory.THistSvc()
+ # thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
+ # acc.addService(thistSvc)
return acc
# with output defaults
diff --git a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
index 814fc31a5b82c9705472a283ba8abfd074ed067b..c8a7df00173d6503683ac47148baefd8bf54da4c 100644
--- a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
+++ b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
@@ -122,7 +122,7 @@ FaserSCT_SpacePoint* TrackerSpacePointMakerTool::makeSCT_SpacePoint(const Tracke
{
if (fabs(lambda1) > 1 + m_stripLengthTolerance)
{
- ATH_MSG_WARNING("Intersection lies outside the bounds of both strips");
+ ATH_MSG_DEBUG("Intersection lies outside the bounds of both strips");
ok = false;
}
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index f990be9458e0ee2ac0a5bdef8b670703f1424c26..14616c6f42dd0e60df67c5b73358da599093130c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -50,6 +50,7 @@ StatusCode CKF2::initialize() {
ATH_CHECK(m_kalmanFitterTool1.retrieve());
ATH_CHECK(m_createTrkTrackTool.retrieve());
ATH_CHECK(m_trackCollection.initialize());
+ // ATH_CHECK(m_allTrackCollection.initialize());
ATH_CHECK(m_eventInfoKey.initialize());
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool.retrieve());
@@ -85,6 +86,9 @@ StatusCode CKF2::execute() {
SG::WriteHandle trackContainer{m_trackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+ // SG::WriteHandle allTrackContainer{m_allTrackCollection, ctx};
+ // std::unique_ptr<TrackCollection> outputAllTracks = std::make_unique<TrackCollection>();
+
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
@@ -157,10 +161,19 @@ StatusCode CKF2::execute() {
// 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");
- }
+ result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps)
+ {
+ try
+ {
+ if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error))
+ {
+ ATH_MSG_WARNING ("Cannot set error state for SCT.");
+ }
+ }
+ catch (...)
+ {
+ ATH_MSG_DEBUG ("SCT error state is locked.");
+ }
}
continue;
}
@@ -183,6 +196,33 @@ StatusCode CKF2::execute() {
// select all tracks with at least 13 heats and with 6 or less shared hits, starting from the best track
// TODO use Gaudi parameters for the number of hits and shared hits
// TODO allow shared hits only in the first station?
+ // std::vector<FaserActsRecMultiTrajectory> rawTrajectories {};
+ // for (auto raw : allTrajectories)
+ // {
+ // rawTrajectories.push_back(raw.trajectory);
+ // }
+ // for (const FaserActsRecMultiTrajectory &traj : rawTrajectories) {
+ // const auto params = traj.trackParameters(traj.tips().front());
+ // ATH_MSG_DEBUG("Fitted parameters (raw)");
+ // ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ // ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ // ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ // ATH_MSG_DEBUG(" charge: " << params.charge());
+ // std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj);
+ // if (track != nullptr) {
+ // m_numberOfSelectedTracks++;
+ // std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, gctx, track.get(), Acts::BoundVector::Zero(), m_isMC, origin);
+ // if (track2) {
+ // outputAllTracks->push_back(std::move(track2));
+ // } else {
+ // outputAllTracks->push_back(std::move(track));
+ // ATH_MSG_WARNING("Re-Fit failed.");
+ // }
+ // } else {
+ // ATH_MSG_WARNING("CKF failed.");
+ // }
+ // }
+
std::vector<FaserActsRecMultiTrajectory> selectedTrajectories {};
while (not allTrajectories.empty()) {
TrajectoryInfo selected = allTrajectories.front();
@@ -224,6 +264,7 @@ StatusCode CKF2::execute() {
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(gctx, selectedTrajectories));
}
+ // ATH_CHECK(allTrackContainer.record(std::move(outputAllTracks)));
ATH_CHECK(trackContainer.record(std::move(outputTracks)));
return StatusCode::SUCCESS;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 170f0353904bef82095146c50bf98d5721ef2138..cd78bf72531931fc3ef7ef3b1b67edc854e46596 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -149,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::WriteHandleKey<TrackCollection> m_allTrackCollection { this, "AllTrackCollection", "CKFAllTrackCollection", "Output all track collection name" };
SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index 952f75287f594e5df31feb2b6a493bfe235a4086..bbae70e3b10be38ff5de9675796bbe509c8f2483 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -130,10 +130,16 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
std::list<Seed> allSeeds;
for (const Seed &seed : seeds) allSeeds.push_back(seed);
+ // allSeeds.sort([](const Seed &left, const Seed &right) {
+ // if (left.size > right.size) return true;
+ // if (left.size < right.size) return false;
+ // if (left.chi2 < right.chi2) return true;
+ // else return false;
+ // });
allSeeds.sort([](const Seed &left, const Seed &right) {
- if (left.size > right.size) return true;
- if (left.size < right.size) return false;
- if (left.chi2 < right.chi2) return true;
+ if (left.stations > right.stations) return true;
+ if (left.stations < right.stations) return false;
+ if (left.chi2/std::max(left.positions.size() + left.fakePositions.size() - left.constraints, 1UL) < right.chi2/std::max(right.positions.size() + right.fakePositions.size() - right.constraints, 1UL)) return true;
else return false;
});
@@ -141,6 +147,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
while (not allSeeds.empty()) {
Seed selected = allSeeds.front();
selectedSeeds.push_back(selected);
+ // allSeeds.remove_if([&](const Seed &p) {
+ // return ((p.clusterSet & selected.clusterSet).count() == p.clusterSet.count());
+ // });
allSeeds.remove_if([&](const Seed &p) {
return (p.size < 10) || ((p.clusterSet & selected.clusterSet).count() > 6);
});
@@ -161,7 +170,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
double origin = !selectedSeeds.empty() ? minSeed->minZ - 10 : 0;
m_targetZPosition = origin;
std::vector<Acts::CurvilinearTrackParameters> initParams {};
+ ATH_MSG_DEBUG("Sorted seed properties:");
for (const Seed &seed : selectedSeeds) {
+ ATH_MSG_DEBUG("seed size: " << seed.size << ", chi2: " << seed.chi2);
initParams.push_back(seed.get_params(origin, cov));
}
@@ -264,13 +275,16 @@ CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
if (segments.size() > 1) {
fakeFit();
+ constraints = 5;
} else {
momentum = 9999999.;
charge = 1;
+ constraints = 2;
}
getChi2();
size = clusters.size();
+ stations = segments.size();
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
index 0467c92620eebd955f0a78f11cfdea3fd0cb0b1d..a4d53c629c24f53e1dfe4488ca7524755f117ad3 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
@@ -66,7 +66,7 @@ private:
double c0, c1, cx, cy, r, chi2, momentum, charge, minZ;
Acts::Vector3 direction;
- size_t size;
+ size_t size, stations, constraints;
Acts::CurvilinearTrackParameters get_params(double origin, Acts::BoundSymMatrix cov) const;
private:
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
index 4928fc917e2b0bfe73e38848939690680fc3c19a..994689010203a269946c3c88392004a61a42adf4 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
@@ -61,9 +61,17 @@ CreateTrkTrackTool::createTrack(const Acts::GeometryContext &gctx, const FaserAc
};
}
}
- const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(clusterOnTrack, parm);
- if (perState) {
- finalTrajectory->insert(finalTrajectory->begin(), perState);
+ double nDoF = state.calibratedSize();
+ const Trk::FitQualityOnSurface* quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
+ const Trk::TrackStateOnSurface* perState = new Trk::TrackStateOnSurface(
+ clusterOnTrack,
+ parm,
+ quality,
+ nullptr,
+ typePattern);
+ if (perState)
+ {
+ finalTrajectory->insert(finalTrajectory->begin(), perState);
}
}
return;