diff --git a/CMakeLists.txt b/CMakeLists.txt
index 33330dc1b94c741e87fea6264e2389cc0293e3d9..7a7399593814cb73cd9d357b18521ae2fbfc87d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -42,6 +42,7 @@ lhcb_add_subdirectories(
     CaloFuture/CaloFutureReco
     CaloFuture/CaloFutureTools
     FT/FTMonitors
+    FT/FTAnalysis
     MicroDST/MicroDSTAlgorithm
     Muon/MuonInterfaces
     Muon/MuonTrackRec
diff --git a/FT/FTAnalysis/CMakeLists.txt b/FT/FTAnalysis/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1679eded1846f3403d0661883d9eeac993186972
--- /dev/null
+++ b/FT/FTAnalysis/CMakeLists.txt
@@ -0,0 +1,21 @@
+###############################################################################
+# (c) Copyright 2023 CERN for the benefit of the LHCb Collaboration           #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+gaudi_add_module(FTAnalysis
+    SOURCES
+        src/FTBeamTimeAnalysis.cpp
+    LINK
+        Gaudi::GaudiAlgLib
+        Gaudi::GaudiCommonSvcLib
+        Gaudi::GaudiKernel
+        LHCb::LHCbKernel
+        LHCb::DAQEventLib
+        LHCb::FTEvent
+)
diff --git a/FT/FTAnalysis/src/FTBeamTimeAnalysis.cpp b/FT/FTAnalysis/src/FTBeamTimeAnalysis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d26e0e29a43be930d7b91324760e190bddd2eb7
--- /dev/null
+++ b/FT/FTAnalysis/src/FTBeamTimeAnalysis.cpp
@@ -0,0 +1,457 @@
+/*****************************************************************************\
+* (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the GNU General Public      *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+*                                                                             *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization  *
+* or submit itself to any jurisdiction.                                       *
+\*****************************************************************************/
+
+#include "Detector/FT/FTChannelID.h"
+#include "Detector/FT/FTConstants.h"
+#include "Event/FTLiteCluster.h"
+#include "Event/ODIN.h"
+#include "Event/TAEUtils.h"
+#include "Gaudi/Accumulators/StaticHistogram.h"
+#include "GaudiAlg/GaudiHistoAlg.h"
+#include "LHCbAlgs/MergingTransformer.h"
+
+#include <map>
+#include <vector>
+
+namespace FTConstants = LHCb::Detector::FT;
+using FTChannelID     = LHCb::Detector::FTChannelID;
+using FTLiteClusters  = LHCb::FTLiteCluster::FTLiteClusters;
+using ODINVector      = Gaudi::Functional::vector_of_const_<LHCb::ODIN const*>;
+using InputVector     = Gaudi::Functional::vector_of_const_<FTLiteClusters const*>;
+
+class FTBeamTimeAnalysis final
+    : public LHCb::Algorithm::MergingConsumer<void( ODINVector const&, InputVector const& )> {
+
+public:
+  FTBeamTimeAnalysis( const std::string& name, ISvcLocator* pSvcLocator );
+
+  StatusCode initialize() override;
+
+  void operator()( ODINVector const&, InputVector const& ) const override;
+
+private:
+  LHCb::TAE::Handler m_taeHandler{ this };
+
+  Gaudi::Property<bool> m_enableChannelHist{ this, "EnableChannelHist", false, "Enable histograms at channel level" };
+
+  Gaudi::Property<bool> m_autoBX{ this, "AutoBX", true,
+                                  "Automatically determine isolated leading or trailing bunch crossing id" };
+
+  Gaudi::Property<int> m_maxClustersPerBX{ this, "MaxClustersPerBX", -1,
+                                           "Maximum of the number of clusters per bunch crossing" };
+  Gaudi::Property<int> m_minClustersPerBX{ this, "MinClustersPerBX", -1,
+                                           "Minimum of the number of clusters per bunch crossing" };
+
+  Gaudi::Property<std::string> m_extraInfo{ this, "ExtraInfo", "", "Extra information shown on the plot titles" };
+
+  Gaudi::Property<unsigned int> m_minBXID{ this, "MinBXID", 0, "Minimum bunch crossing ID" };
+  Gaudi::Property<unsigned int> m_maxBXID{ this, "MaxBXID", 3563, "Maximum bunch crossing ID" };
+
+  Gaudi::Property<unsigned int> m_minStep{ this, "MinStep", 0, "Minimum step number" };
+  Gaudi::Property<unsigned int> m_numStep{ this, "NumStep", 1, "Number of steps" };
+
+  Gaudi::Property<unsigned int> m_taeHalfWindow{ this, "TaeHalfWindow", 7, "TAE half window" };
+
+  // Manually determine isolated or leading or trailing
+  // To be used when AutoBX is disabled
+  Gaudi::Property<std::vector<unsigned int>> m_isobxid{ this, "IsolatedBXID", {}, "Isolated BXID" };
+  Gaudi::Property<std::vector<unsigned int>> m_leadbxid{ this, "LeadingBXID", {}, "Leading BXID" };
+  Gaudi::Property<std::vector<unsigned int>> m_trailbxid{ this, "TrailingBXID", {}, "Trailing BXID" };
+
+  using Histogram1D = Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float>;
+  using Histogram2D = Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float>;
+  using Axis        = Gaudi::Accumulators::Axis<float>;
+  using H2DArg      = Gaudi::Accumulators::HistoInputType<float, 2>;
+
+  struct TaeHists {
+    TaeHists( const FTBeamTimeAnalysis* parent, const std::string& nameprefix, const std::string namesuffix,
+              const std::string& extrainfo, bool enable_channelhist, unsigned int tae_half_window )
+        : m_taehist_perlink{
+              parent, fmt::format( "{}TAE/LinkLevel/{}", nameprefix, namesuffix ),
+              fmt::format( "{} TAE {} {};Datalink Index;Relative BXID;Number of clusters", nameprefix, namesuffix,
+                           extrainfo ),
+              std::tuple{ Axis{ FTConstants::nSiPMsTotal, -0.5f, FTConstants::nSiPMsTotal - 0.5f },
+                          Axis{ 2 * tae_half_window + 1, -1.0f * static_cast<float>( tae_half_window ) - 0.5f,
+                                static_cast<float>( tae_half_window ) + 0.5f } } } {
+      if ( enable_channelhist ) {
+        for ( unsigned int i_station : { 1, 2, 3 } ) {
+          for ( unsigned int i_layer : { 0, 1, 2, 3 } ) {
+            for ( unsigned int i_quarter : { 0, 1, 2, 3 } ) {
+              auto tmp_channel_id = FTChannelID( FTChannelID::StationID( i_station ), FTChannelID::LayerID( i_layer ),
+                                                 FTChannelID::QuarterID( i_quarter ), FTChannelID::ModuleID( 0 ),
+                                                 FTChannelID::MatID( 0 ), 0, 0 );
+              unsigned int global_quarter_idx = tmp_channel_id.globalQuarterIdx();
+
+              m_taehist_perchannel.emplace(
+                  std::piecewise_construct, std::forward_as_tuple( global_quarter_idx ),
+                  std::forward_as_tuple(
+                      parent,
+                      fmt::format( "{}TAE/ChannelLevel/T{}L{}Q{}/{}", nameprefix, i_station, i_layer, i_quarter,
+                                   namesuffix ),
+                      fmt::format( "{} TAE {} {};Channel index;Relative BXID;Number of clusters", nameprefix,
+                                   namesuffix, extrainfo ),
+                      std::tuple{ Axis{ FTConstants::nMaxChannelsPerQuarter, -0.5f,
+                                        FTConstants::nMaxChannelsPerQuarter - 0.5f },
+                                  Axis{ 2 * tae_half_window + 1, -1.0f * static_cast<float>( tae_half_window ) - 0.5f,
+                                        static_cast<float>( tae_half_window ) + 0.5f } } ) );
+            }
+          }
+        }
+      }
+    }
+
+    void FillCluster( int relative_bxid, const FTChannelID& chan_id ) {
+      auto global_sipm_idx = chan_id.globalSipmIdx();
+      m_taehist_perlink[H2DArg{ global_sipm_idx, relative_bxid }] += 1;
+      if ( !m_taehist_perchannel.empty() ) {
+        auto global_quarter_idx           = chan_id.globalQuarterIdx();
+        auto local_channel_idx_in_quarter = chan_id.localChannelIdx_quarter();
+        m_taehist_perchannel.at( global_quarter_idx )[H2DArg{ local_channel_idx_in_quarter, relative_bxid }] += 1;
+      }
+    }
+
+    Histogram2D m_taehist_perlink;
+    // global quarter index -> histogram
+    std::map<unsigned int, Histogram2D> m_taehist_perchannel;
+  };
+
+  struct StepTaeHistGroup {
+    StepTaeHistGroup( const FTBeamTimeAnalysis* parent, const std::string& nameprefix, const std::string& extrainfo,
+                      bool enable_channelhist, unsigned int min_step_index, unsigned int max_step_index,
+                      unsigned tae_half_window ) {
+      for ( unsigned int i_step = min_step_index; i_step <= max_step_index; i_step++ ) {
+        m_steptaehists.emplace( std::piecewise_construct, std::forward_as_tuple( i_step ),
+                                std::forward_as_tuple( parent, nameprefix, fmt::format( "Step{}", i_step ), extrainfo,
+                                                       enable_channelhist, tae_half_window ) );
+      }
+    }
+
+    void FillCluster( unsigned int step_index, int relative_bxid, const FTChannelID& chan_id ) {
+      m_steptaehists.at( step_index ).FillCluster( relative_bxid, chan_id );
+    }
+
+    std::map<unsigned int, TaeHists> m_steptaehists;
+  };
+
+  enum class TAEType { Others, IsoBX, LeadBX, TrailBX };
+
+  // information concerning all bxid, including TAE and nonTAE
+  mutable std::optional<Histogram1D> m_hist_evtperbx;
+  mutable std::optional<Histogram1D> m_hist_clusterperbx;
+  mutable std::optional<Histogram2D> m_hist_singleevtclusterperbx;
+
+  // information about events which failed the number of clusters per event cut
+  mutable std::optional<Histogram1D> m_hist_highocc_evt_perbx;
+  mutable std::optional<Histogram1D> m_hist_lowocc_evt_perbx;
+
+  // summary information, for IsoBX, LeadBX, TrailBX TAE combined
+  mutable std::optional<Histogram1D> m_hist_clusterperlink;
+
+  // nonstep TAE histogram
+  mutable std::map<TAEType, TaeHists> m_nonstep_tae_hists;
+
+  // Only for TAE and step run
+  mutable std::map<TAEType, Histogram1D> m_hist_evtperstep;
+  mutable std::map<TAEType, Histogram1D> m_hist_clusterperstep;
+
+  // step TAE histograms
+  mutable std::map<TAEType, StepTaeHistGroup> m_step_tae_hists;
+
+  // Determining bxid flavor based on EventType field in ODIN bank
+  bool AroundIsolatedBX( std::uint16_t event_type ) const {
+    // Isolated bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_08 ) ) { return true; }
+    // Empty before isolated bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_11 ) ) { return true; }
+    // Empty after isolated bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_12 ) ) { return true; }
+    return false;
+  }
+
+  // Determining bxid flavor based on EventType field in ODIN bank
+  bool AroundLeadingBX( std::uint16_t event_type ) const {
+    // Leading bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_09 ) ) { return true; }
+    // Empty before Leading bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_13 ) ) { return true; }
+    return false;
+  }
+
+  // Determining bxid flavor based on EventType field in ODIN bank
+  bool AroundTrailingBX( std::uint16_t event_type ) const {
+    // Trailing bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_10 ) ) { return true; }
+    // Empty afer Trailing bxid
+    if ( event_type & static_cast<std::uint16_t>( LHCb::ODIN::EventTypes::et_bit_14 ) ) { return true; }
+    return false;
+  }
+
+  TAEType AutoDetermineTaeType( const LHCb::ODIN& odin ) const {
+    std::uint16_t event_type = odin.eventType();
+    if ( AroundIsolatedBX( event_type ) ) { return TAEType::IsoBX; }
+    if ( AroundLeadingBX( event_type ) ) { return TAEType::LeadBX; }
+    if ( AroundTrailingBX( event_type ) ) { return TAEType::TrailBX; }
+    return TAEType::Others;
+  }
+
+  TAEType ManualDetermineTaeType( int evt_bxid, int relative_bxid ) const {
+    if ( !m_isobxid.empty() ) {
+      for ( int test_center_bxid : m_isobxid.value() ) {
+        int test_target_bxid = test_center_bxid + relative_bxid;
+        if ( evt_bxid == test_target_bxid ) { return TAEType::IsoBX; }
+      }
+    }
+    if ( !m_leadbxid.empty() ) {
+      for ( int test_center_bxid : m_leadbxid.value() ) {
+        int test_target_bxid = test_center_bxid + relative_bxid;
+        if ( evt_bxid == test_target_bxid ) { return TAEType::LeadBX; }
+      }
+    }
+    if ( !m_trailbxid.empty() ) {
+      for ( int test_center_bxid : m_trailbxid.value() ) {
+        int test_target_bxid = test_center_bxid + relative_bxid;
+        if ( evt_bxid == test_target_bxid ) { return TAEType::TrailBX; }
+      }
+    }
+    return TAEType::Others;
+  }
+};
+
+DECLARE_COMPONENT( FTBeamTimeAnalysis )
+
+FTBeamTimeAnalysis::FTBeamTimeAnalysis( const std::string& name, ISvcLocator* pSvcLocator )
+    : LHCb::Algorithm::MergingConsumer<void( ODINVector const&, InputVector const& )>(
+          name, pSvcLocator, { KeyValues{ "ODINVector", {} }, KeyValues{ "InputVector", {} } } ){};
+
+StatusCode FTBeamTimeAnalysis::initialize() {
+  return LHCb::Algorithm::MergingConsumer<void( ODINVector const&, InputVector const& )>::initialize().andThen( [&] {
+    m_hist_evtperbx.emplace(
+        this, "EvtPerBX", fmt::format( "Number of events per BX {};BXID;Number of events", m_extraInfo.value() ),
+        Axis{ m_maxBXID.value() - m_minBXID.value() + 1, m_minBXID.value() - 0.5f, m_maxBXID.value() + 0.5f } );
+    m_hist_clusterperbx.emplace(
+        this, "ClusterPerBX",
+        fmt::format( "Number of clusters per BX {};BXID;Number of clusters", m_extraInfo.value() ),
+        Axis{ m_maxBXID.value() - m_minBXID.value() + 1, m_minBXID.value() - 0.5f, m_maxBXID.value() + 0.5f } );
+    m_hist_singleevtclusterperbx.emplace(
+        this, "SingleEvtClusterPerBX",
+        fmt::format(
+            "Singe event number of clusters per BX {};BXID;Single event number of clusters;Number of events per bin",
+            m_extraInfo.value() ),
+        std::tuple{
+            Axis{ m_maxBXID.value() - m_minBXID.value() + 1, m_minBXID.value() - 0.5f, m_maxBXID.value() + 0.5f },
+            Axis{ 101, -50.0, 10050.0 } } );
+
+    m_hist_clusterperlink.emplace(
+        this, "TAENumClusterPerLink",
+        fmt::format( "Number of clusters per link {};Datalink Index;Number of clusters", m_extraInfo.value() ),
+        Axis{ FTConstants::nSiPMsTotal, -0.5f, FTConstants::nSiPMsTotal - 0.5f } );
+
+    if ( m_maxClustersPerBX.value() > 0 ) {
+      m_hist_highocc_evt_perbx.emplace(
+          this, "HighOccupancyEvts",
+          fmt::format( "Number of high occupancy events per BX {};BXID;Number of events", m_extraInfo.value() ),
+          Axis{ m_maxBXID.value() - m_minBXID.value() + 1, m_minBXID.value() - 0.5f, m_maxBXID.value() + 0.5f } );
+    }
+    if ( m_minClustersPerBX.value() > 0 ) {
+      m_hist_lowocc_evt_perbx.emplace(
+          this, "LowOccupancyEvts",
+          fmt::format( "Number of low occupancy events per BX {};BXID;Number of events", m_extraInfo.value() ),
+          Axis{ m_maxBXID.value() - m_minBXID.value() + 1, m_minBXID.value() - 0.5f, m_maxBXID.value() + 0.5f } );
+    }
+
+    bool fill_isobxid   = ( m_autoBX.value() ) || ( !m_isobxid.empty() );
+    bool fill_leadbxid  = ( m_autoBX.value() ) || ( !m_leadbxid.empty() );
+    bool fill_trailbxid = ( m_autoBX.value() ) || ( !m_trailbxid.empty() );
+
+    if ( 0 == m_numStep.value() ) {
+      if ( fill_isobxid ) {
+        m_nonstep_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::IsoBX ),
+                                     std::forward_as_tuple( this, "Isobx", "Nonstep", m_extraInfo.value(),
+                                                            m_enableChannelHist.value(), m_taeHalfWindow.value() ) );
+      }
+      if ( fill_leadbxid ) {
+        m_nonstep_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::LeadBX ),
+                                     std::forward_as_tuple( this, "Leadbx", "Nonstep", m_extraInfo.value(),
+                                                            m_enableChannelHist.value(), m_taeHalfWindow.value() ) );
+      }
+      if ( fill_trailbxid ) {
+        m_nonstep_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::TrailBX ),
+                                     std::forward_as_tuple( this, "Trailbx", "Nonstep", m_extraInfo.value(),
+                                                            m_enableChannelHist.value(), m_taeHalfWindow.value() ) );
+      }
+    } else {
+      unsigned int max_step = m_minStep.value() + m_numStep.value() - 1;
+      if ( fill_isobxid ) {
+        m_hist_evtperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::IsoBX ),
+            std::forward_as_tuple(
+                this, "IsobxNumEvtPerStep",
+                fmt::format( "Number of events (Isolated BX) per step {};Step;Number of events", m_extraInfo.value() ),
+                Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_hist_clusterperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::IsoBX ),
+            std::forward_as_tuple( this, "IsobxNumClusterPerStep",
+                                   fmt::format( "Number of clusters (Isolated BX) per step {};Step;Number of clusters",
+                                                m_extraInfo.value() ),
+                                   Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_step_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::IsoBX ),
+                                  std::forward_as_tuple( this, "Isobx", m_extraInfo.value(),
+                                                         m_enableChannelHist.value(), m_minStep.value(), max_step,
+                                                         m_taeHalfWindow.value() ) );
+      }
+      if ( fill_leadbxid ) {
+        m_hist_evtperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::LeadBX ),
+            std::forward_as_tuple(
+                this, "LeadbxNumEvtPerStep",
+                fmt::format( "Number of events (Leading BX) per step {};Step;Number of events", m_extraInfo.value() ),
+                Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_hist_clusterperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::LeadBX ),
+            std::forward_as_tuple( this, "LeadbxNumClusterPerStep",
+                                   fmt::format( "Number of clusters (Leading BX) per step {};Step;Number of clusters",
+                                                m_extraInfo.value() ),
+                                   Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_step_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::LeadBX ),
+                                  std::forward_as_tuple( this, "Leadbx", m_extraInfo.value(),
+                                                         m_enableChannelHist.value(), m_minStep.value(), max_step,
+                                                         m_taeHalfWindow.value() ) );
+      }
+      if ( fill_trailbxid ) {
+        m_hist_evtperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::TrailBX ),
+            std::forward_as_tuple(
+                this, "TrailbxNumEvtPerStep",
+                fmt::format( "Number of events (Trailing BX) per step {};Step;Number of events", m_extraInfo.value() ),
+                Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_hist_clusterperstep.emplace(
+            std::piecewise_construct, std::forward_as_tuple( TAEType::TrailBX ),
+            std::forward_as_tuple( this, "TrailbxNumClusterPerStep",
+                                   fmt::format( "Number of clusters (Trailing BX) per step {};Step;Number of clusters",
+                                                m_extraInfo.value() ),
+                                   Axis{ m_numStep.value(), m_minStep.value() - 0.5f, max_step + 0.5f } ) );
+        m_step_tae_hists.emplace( std::piecewise_construct, std::forward_as_tuple( TAEType::TrailBX ),
+                                  std::forward_as_tuple( this, "Trailbx", m_extraInfo.value(),
+                                                         m_enableChannelHist.value(), m_minStep.value(), max_step,
+                                                         m_taeHalfWindow.value() ) );
+      }
+    }
+
+    return StatusCode::SUCCESS;
+  } );
+}
+
+void FTBeamTimeAnalysis::operator()( ODINVector const& odinVector, InputVector const& inputVector ) const {
+
+  // without specify the taeHandler nSide, try to get all TAE events
+  auto taeEvents = m_taeHandler.arrangeTAE( odinVector, inputVector );
+  if ( taeEvents.empty() ) { return; }
+
+  // If there is requirement on the number of clusters per event
+  // First find the maximum number of clusters in the several consecutive events in TAE
+  // Then keep or discard the whole TAE according to that number
+  // Also determining the whole TAE array as isolated, leading or trailing
+  int     max_occ_per_tae_array = 0;
+  TAEType arr_tae_type          = TAEType::Others;
+  bool    found_iso_in_arr      = false;
+  bool    found_lead_in_arr     = false;
+  bool    found_trail_in_arr    = false;
+  for ( const auto& element : taeEvents ) {
+    int                   relative_bxid       = element.first;
+    LHCb::ODIN const&     odin                = element.second.first;
+    FTLiteClusters const& clusters            = element.second.second;
+    int                   bx_id               = odin.bunchId();
+    int                   num_clusters_per_bx = static_cast<int>( clusters.size() );
+    if ( num_clusters_per_bx > max_occ_per_tae_array ) { max_occ_per_tae_array = num_clusters_per_bx; }
+    // The histogram single event number of clusters is filled without considering the number of clusters per event
+    // limit
+    ( *m_hist_singleevtclusterperbx )[H2DArg{ bx_id, num_clusters_per_bx }] += 1;
+
+    // Determine the TAEType of the single event
+    TAEType evt_tae_type = TAEType::Others;
+    if ( m_autoBX.value() ) {
+      evt_tae_type = AutoDetermineTaeType( odin );
+    } else {
+      evt_tae_type = ManualDetermineTaeType( bx_id, relative_bxid );
+    }
+    if ( TAEType::IsoBX == evt_tae_type ) { found_iso_in_arr = true; }
+    if ( TAEType::LeadBX == evt_tae_type ) { found_lead_in_arr = true; }
+    if ( TAEType::TrailBX == evt_tae_type ) { found_trail_in_arr = true; }
+  }
+
+  // If none of the Iso, Lead, Trail were found or more than one of them were found, determine the TAE array as "Others"
+  if ( found_iso_in_arr && ( !found_lead_in_arr ) && ( !found_trail_in_arr ) ) {
+    arr_tae_type = TAEType::IsoBX;
+  } else if ( ( !found_iso_in_arr ) && found_lead_in_arr && ( !found_trail_in_arr ) ) {
+    arr_tae_type = TAEType::LeadBX;
+  } else if ( ( !found_iso_in_arr ) && ( !found_lead_in_arr ) && found_trail_in_arr ) {
+    arr_tae_type = TAEType::TrailBX;
+  } else {
+    arr_tae_type = TAEType::Others;
+  }
+
+  int limit_max_clusters_per_bx = m_maxClustersPerBX.value();
+  int limit_min_clusters_per_bx = m_minClustersPerBX.value();
+
+  bool evt_occ_too_high = false;
+  if ( limit_max_clusters_per_bx >= 0 ) {
+    if ( limit_max_clusters_per_bx < max_occ_per_tae_array ) { evt_occ_too_high = true; }
+  }
+  bool evt_occ_too_low = false;
+  if ( limit_min_clusters_per_bx >= 0 ) {
+    if ( limit_min_clusters_per_bx > max_occ_per_tae_array ) { evt_occ_too_low = true; }
+  }
+  // Cut away those beyond the number of clusters per event limit
+  if ( evt_occ_too_high || evt_occ_too_low ) {
+    for ( const auto& element : taeEvents ) {
+      LHCb::ODIN const&     odin                = element.second.first;
+      FTLiteClusters const& clusters            = element.second.second;
+      int                   bx_id               = odin.bunchId();
+      int                   num_clusters_per_bx = static_cast<int>( clusters.size() );
+      if ( limit_max_clusters_per_bx >= 0 ) {
+        if ( limit_max_clusters_per_bx < num_clusters_per_bx ) { ( *m_hist_highocc_evt_perbx )[bx_id] += 1; }
+      }
+      if ( limit_min_clusters_per_bx >= 0 ) {
+        if ( limit_min_clusters_per_bx > num_clusters_per_bx ) { ( *m_hist_lowocc_evt_perbx )[bx_id] += 1; }
+      }
+    }
+    return;
+  }
+
+  for ( const auto& element : taeEvents ) {
+    int                   relative_bxid = element.first;
+    LHCb::ODIN const&     odin          = element.second.first;
+    int                   bx_id         = odin.bunchId();
+    FTLiteClusters const& clusters      = element.second.second;
+
+    ( *m_hist_evtperbx )[bx_id] += 1;
+    ( *m_hist_clusterperbx )[bx_id] += clusters.size();
+
+    if ( TAEType::Others != arr_tae_type ) {
+      if ( 0 == m_numStep.value() ) {
+        for ( const auto& cluster : clusters.range() ) {
+          FTChannelID chan_id = cluster.channelID();
+          ( *m_hist_clusterperlink )[chan_id.globalSipmIdx()] += 1;
+          m_nonstep_tae_hists.at( arr_tae_type ).FillCluster( relative_bxid, chan_id );
+        }
+      } else {
+        auto step_idx = odin.calibrationStep();
+        m_hist_evtperstep.at( arr_tae_type )[step_idx] += 1;
+        m_hist_clusterperstep.at( arr_tae_type )[step_idx] += clusters.size();
+        for ( const auto& cluster : clusters.range() ) {
+          FTChannelID chan_id = cluster.channelID();
+          ( *m_hist_clusterperlink )[chan_id.globalSipmIdx()] += 1;
+          m_step_tae_hists.at( arr_tae_type ).FillCluster( step_idx, relative_bxid, chan_id );
+        }
+      }
+    }
+  }
+}