From e18e8d3fee07563ae97a1132ce13ba2b0ac75244 Mon Sep 17 00:00:00 2001
From: FASER Reco <faserrec@lxplus780.cern.ch>
Date: Sat, 22 Apr 2023 01:09:48 +0200
Subject: [PATCH] Fixes merged from p0010

---
 .../Reconstruction/scripts/faser_reco.py      | 48 ++++++++-----
 LHCData/LHCDataAlgs/src/LHCDataAlg.cxx        | 69 ++++++++++++++++---
 2 files changed, 88 insertions(+), 29 deletions(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 21e1d51a..744aafa0 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/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx b/LHCData/LHCDataAlgs/src/LHCDataAlg.cxx
index de226985..13229039 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;
 
 }
-- 
GitLab