diff --git a/Dumpers/BinaryDumpers/options/dump_geometry.py b/Dumpers/BinaryDumpers/options/dump_geometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..005f42edb8d7f27d1f9f698e5d63802bdc7125ca
--- /dev/null
+++ b/Dumpers/BinaryDumpers/options/dump_geometry.py
@@ -0,0 +1,103 @@
+###############################################################################
+# (c) Copyright 2000-2018 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.                                       #
+###############################################################################
+from GaudiConf import IOHelper
+from GaudiPython.Bindings import AppMgr, gbl
+from Configurables import LHCbApp, CondDB
+from Configurables import GaudiSequencer
+from Configurables import FTRawBankDecoder
+from Configurables import (PrStoreUTHit, PrStoreFTHit)
+from Configurables import ApplicationMgr
+from Configurables import HistogramPersistencySvc
+from Configurables import (AuditorSvc, SequencerTimerTool)
+from Configurables import (DumpUTGeometry, DumpFTGeometry, DumpMuonTable,
+                           DumpUTLookupTables, DumpMuonGeometry,
+                           DumpCaloGeometry, DumpVPGeometry, DumpMagneticField,
+                           DumpBeamline)
+from Configurables import RootHistCnv__PersSvc
+from Configurables import IODataManager
+from Configurables import createODIN
+from Configurables import TestMuonTable
+
+app = LHCbApp(
+    DataType="Upgrade",
+    EvtMax=50,
+    Simulation=True,
+    DDDBtag="dddb-20171122",
+    CondDBtag="sim-20180530-vc-md100")
+
+# Upgrade DBs
+CondDB().Upgrade = True
+
+dec_seq = GaudiSequencer("DecodingSeq")
+dec_seq.Members = [createODIN()]
+
+ecal_location = "/dd/Structure/LHCb/DownstreamRegion/Ecal"
+hcal_location = "/dd/Structure/LHCb/DownstreamRegion/Hcal"
+
+ecal_geom = DumpCaloGeometry("DumpEcal", Location=ecal_location)
+hcal_geom = DumpCaloGeometry("DumpHcal", Location=hcal_location)
+
+# Add the service that will dump the UT and FT geometry
+ApplicationMgr().ExtSvc += [
+    DumpMagneticField(),
+    DumpBeamline(),
+    DumpVPGeometry(),
+    DumpUTGeometry(),
+    DumpFTGeometry(),
+    DumpMuonGeometry(),
+    DumpMuonTable(),
+    DumpUTLookupTables(), ecal_geom, hcal_geom
+]
+
+ApplicationMgr().TopAlg = [dec_seq]
+
+# Some extra stuff for timing table
+ApplicationMgr().ExtSvc += ['ToolSvc', 'AuditorSvc']
+ApplicationMgr().AuditAlgorithms = True
+AuditorSvc().Auditors += ['TimingAuditor']
+SequencerTimerTool().OutputLevel = 4
+
+# Some extra stuff to save histograms
+ApplicationMgr().HistogramPersistency = "ROOT"
+RootHistCnv__PersSvc('RootHistCnv').ForceAlphaIds = True
+HistogramPersistencySvc().OutputFile = "dump-histos.root"
+
+# No error messages when reading MDF
+IODataManager().DisablePFNWarning = True
+
+IOHelper('MDF').inputFiles(
+    [
+        # SciFi v5, minbias
+        'PFN:/data/bfys/raaij/upgrade_scifi_v5_uncompressed/upgrade_mc_minbias_scifi_v5_%03d.mdf'
+        % i for i in range(5)
+    ],
+    clear=True)
+
+# GaudiPython
+gaudi = AppMgr()
+gaudi.initialize()
+
+TES = gaudi.evtSvc()
+det = gaudi.detSvc()
+ecal = det[ecal_location]
+hcal = det[hcal_location]
+
+ecal_type = gbl.LHCb.RawBank.EcalPacked
+hcal_type = gbl.LHCb.RawBank.HcalPacked
+
+while True:
+    gaudi.run(1)
+    if not TES['/Event']:
+        break
+    raw_event = TES['DAQ/RawEvent']
+    ecal_banks = raw_event.banks(ecal_type)
+    hcal_banks = raw_event.banks(hcal_type)
+    break
diff --git a/Rec/Allen/CMakeLists.txt b/Rec/Allen/CMakeLists.txt
index 4d651bd4080ced818dc2ac7565596301edc5f666..4246396a5b94975fb93970f3a4a93de9e0a0dc51 100644
--- a/Rec/Allen/CMakeLists.txt
+++ b/Rec/Allen/CMakeLists.txt
@@ -57,6 +57,7 @@ include_directories(${CMAKE_SOURCE_DIR}/device/vertex_fit/common/include)
 include_directories(${CMAKE_SOURCE_DIR}/device/vertex_fit/vertex_fitter/include)
 include_directories(${CMAKE_SOURCE_DIR}/device/selections/line_types/include)
 include_directories(${CMAKE_SOURCE_DIR}/device/selections/Hlt1/include)
+include_directories(${CMAKE_SOURCE_DIR}/device/calo/clustering/include)
 include_directories(${CMAKE_SOURCE_DIR}/device/calo/event_model/include)
 include_directories(${CMAKE_SOURCE_DIR}/backend/include)
 include_directories(${CMAKE_SOURCE_DIR}/test/contracts/include)
diff --git a/Rec/Allen/src/AllenCaloToCaloClusters.cpp b/Rec/Allen/src/AllenCaloToCaloClusters.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aa82b7ab92bc45cdd271f4d30b6918d62452f8ee
--- /dev/null
+++ b/Rec/Allen/src/AllenCaloToCaloClusters.cpp
@@ -0,0 +1,189 @@
+/***************************************************************************** \
+ * (c) Copyright 2000-2018 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.                                       *
+\*****************************************************************************/
+
+/**
+ * Convert AllenCalo to CaloCluster v2
+ *
+ * author Dorothea vom Bruch
+ *
+ */
+
+#include "AllenCaloToCaloClusters.h"
+
+DECLARE_COMPONENT(AllenCaloToCaloClusters)
+
+AllenCaloToCaloClusters::AllenCaloToCaloClusters(const std::string& name, ISvcLocator* pSvcLocator) :
+  MultiTransformer(
+    name,
+    pSvcLocator,
+    // Inputs
+    {KeyValue {"AllenOutput", "Allen/Out/HostBuffers"}},
+    // Outputs
+    {KeyValue {"AllenEcalClusters", "Allen/Calo/EcalCluster"},
+     KeyValue {"AllenHcalClusters", "Allen/Calo/HcalCluster"}})
+{}
+
+std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters> AllenCaloToCaloClusters::operator()(
+  const HostBuffers& host_buffers) const
+{
+  // avoid long names
+  // using namespace LHCb::CaloDataFunctor;
+  LHCb::Event::Calo::Clusters EcalClusters;
+  LHCb::Event::Calo::Clusters HcalClusters;
+  // Make the clusters
+  const unsigned i_event = 0;
+  const unsigned number_of_events = 1;
+
+  unsigned number_of_ecal_clusters =
+    host_buffers.host_ecal_cluster_offsets[number_of_events] - host_buffers.host_ecal_cluster_offsets[i_event];
+  unsigned number_of_hcal_clusters =
+    host_buffers.host_hcal_cluster_offsets[number_of_events] - host_buffers.host_hcal_cluster_offsets[i_event];
+  CaloCluster* ecal_clusters = (CaloCluster*) (host_buffers.host_ecal_clusters.data());
+  CaloCluster* hcal_clusters = (CaloCluster*) (host_buffers.host_hcal_clusters.data());
+
+  if (msgLevel(MSG::DEBUG)) {
+    debug() << "Number of Ecal clusters to convert = " << number_of_ecal_clusters << endmsg;
+    debug() << "Number of Hcal clusters to convert = " << number_of_hcal_clusters << endmsg;
+  }
+
+  EcalClusters.reserveForEntries(number_of_ecal_clusters);
+  HcalClusters.reserveForEntries(number_of_hcal_clusters);
+
+  // Loop over Allen Ecal clusters and convert them
+  // Don't need to access them with offset since one event is processed at a time
+  int16_t iFirstEntry = 0;
+  for (unsigned i = 0; i < number_of_ecal_clusters; i++) {
+    const auto& cluster = ecal_clusters[i];
+
+    auto seedCellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.center_id);
+
+    if (msgLevel(MSG::DEBUG) && i < 50 && i % 5 == 0) {
+      if (LHCb::Calo::isValid(seedCellID))
+        debug() << "ECAL CellID " << seedCellID << " corresponding to dense ID " << cluster.center_id << " is invalid!"
+                << endmsg;
+      debug() << " \t ECAL center_id = " << cluster.center_id << " cellID: " << seedCellID << ", e = " << cluster.e
+              << ", x = " << cluster.x << ", y = " << cluster.y;
+      for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j)
+        debug() << " " << cluster.digits[j];
+      debug() << endmsg;
+    }
+
+    // Add the all digits, marking the seed ones
+    EcalClusters.emplace_back(
+      seedCellID,
+      cluster.e,
+      1.0,
+      LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, LHCb::CaloDigitStatus::Mask::SeedCell});
+    for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) {
+      auto cellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.digits[j]);
+      EcalClusters.emplace_back(
+        cellID,
+        0.,
+        1.0,
+        LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy,
+                                       LHCb::CaloDigitStatus::Mask::OwnedCell});
+    }
+    EcalClusters.emplace_back(
+      seedCellID,
+      LHCb::Event::Calo::Clusters::Type::Area3x3,
+      {iFirstEntry, Calo::Constants::max_neighbours + 1},
+      cluster.e,
+      {cluster.x, cluster.y, 0.});
+
+    iFirstEntry += Calo::Constants::max_neighbours + 1; // seed digit+ associated digits making the cluster
+  }
+
+  // Loop over Allen Hcal clusters and convert them
+
+  iFirstEntry = 0;
+  for (unsigned i = 0; i < number_of_hcal_clusters; i++) {
+    const auto& cluster = hcal_clusters[i];
+
+    auto seedCellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.center_id);
+
+    if (msgLevel(MSG::DEBUG) && i < 50 && i % 5 == 0) {
+      if (LHCb::Calo::isValid(seedCellID))
+        debug() << "HCAL CellID " << seedCellID << " corresponding to dense ID " << cluster.center_id << " is invalid!"
+                << endmsg;
+      debug() << " \t HCAL center_id = " << cluster.center_id << " cellID: " << seedCellID << ", e = " << cluster.e
+              << ", x = " << cluster.x << ", y = " << cluster.y;
+      for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j)
+        debug() << " " << cluster.digits[j];
+      debug() << endmsg;
+    }
+
+    // Add the all digits, marking the seed ones
+    HcalClusters.emplace_back(
+      seedCellID,
+      cluster.e,
+      LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy, LHCb::CaloDigitStatus::Mask::SeedCell});
+    for (unsigned j = 0; j < Calo::Constants::max_neighbours; ++j) {
+      auto cellID = LHCb::Calo::DenseIndex::details::toCellID(cluster.digits[j]);
+      HcalClusters.emplace_back(
+        cellID,
+        .0,
+        LHCb::CaloDigitStatus::Status {LHCb::CaloDigitStatus::Mask::UseForEnergy,
+                                       LHCb::CaloDigitStatus::Mask::OwnedCell});
+    }
+    LHCb::Event::Calo::v2::Clusters::which_entries_t entries {iFirstEntry, Calo::Constants::max_neighbours + 1};
+    Gaudi::XYZPoint position {cluster.x, cluster.y, 0.};
+
+    HcalClusters.emplace_back(seedCellID, LHCb::Event::Calo::Clusters::Type::Area3x3, entries, cluster.e, position);
+
+    iFirstEntry += Calo::Constants::max_neighbours + 1; // seed digit+ associated digits making the cluster
+  }
+  if (msgLevel(MSG::DEBUG)) {
+    debug() << "Number of ecal seed clusters: " << EcalClusters.size() << endmsg;
+    uint i = 0;
+    for (const auto& Cluster : EcalClusters) {
+      auto cellID = Cluster.cellID();
+      const double e = Cluster.energy();
+      const double x = Cluster.position().x();
+      const double y = Cluster.position().y();
+      const double z = Cluster.position().z();
+      // const double et = e * sqrt( x * x + y * y ) / sqrt( x * x + y * y + z * z );
+
+      if (i % 5 == 0) {
+        debug() << "Ecal cellID: " << cellID << " energy = " << e << ", x = " << x << ", y = " << y << ", z = " << z
+                << endmsg;
+        auto digits = Cluster.entries();
+        for (const auto& digit : digits)
+          debug() << "     cellID: " << digit.cellID() << " energy: " << digit.energy()
+                  << " fraction: " << digit.fraction() << " Status: " << digit.status() << endmsg;
+      }
+      if (i > 50) break;
+      ++i;
+    }
+    debug() << "Number of hcal seed clusters: " << HcalClusters.size() << endmsg;
+    i = 0;
+    for (const auto& Cluster : HcalClusters) {
+      auto cellID = Cluster.cellID();
+      const double e = Cluster.energy();
+      const double x = Cluster.position().x();
+      const double y = Cluster.position().y();
+      const double z = Cluster.position().z();
+      // const double et = e * sqrt( x * x + y * y ) / sqrt( x * x + y * y + z * z );
+
+      if (i % 5 == 0) {
+        debug() << "Hcal cellID: " << cellID << " energy = " << e << ", x = " << x << ", y = " << y << ", z = " << z
+                << endmsg;
+        auto digits = Cluster.entries();
+        for (const auto& digit : digits)
+          debug() << "     cellID: " << digit.cellID() << " energy: " << digit.energy()
+                  << " fraction: " << digit.fraction() << " Status: " << digit.status() << endmsg;
+      }
+      if (i > 50) break;
+      ++i;
+    }
+  }
+
+  return {EcalClusters, HcalClusters};
+}
diff --git a/Rec/Allen/src/AllenCaloToCaloClusters.h b/Rec/Allen/src/AllenCaloToCaloClusters.h
new file mode 100644
index 0000000000000000000000000000000000000000..5c43b659b05f4759216fdcbf904dca9cc8915b93
--- /dev/null
+++ b/Rec/Allen/src/AllenCaloToCaloClusters.h
@@ -0,0 +1,47 @@
+/***************************************************************************** \
+ * (c) Copyright 2000-2018 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.                                       *
+\*****************************************************************************/
+#ifndef ALLENCALOTOCALOCLUSTERS_H
+#define ALLENCALOTOCALOCLUSTERS_H
+
+// Gaudi
+#include "GaudiAlg/Transformer.h"
+
+// LHCb
+#include "Event/Track.h"
+
+// Allen
+#include "HostBuffers.cuh"
+#include "Logger.h"
+#include "VeloConsolidated.cuh"
+#include "CaloCluster.cuh"
+#include "CaloConstants.cuh"
+//#include "Event/CaloCluster.h"
+#include "Event/CaloHypos_v2.h"
+#include "Event/CaloPosition.h"
+#include "Event/CaloClusters_v2.h"
+#include "Kernel/CaloCellID.h"
+#include "GaudiKernel/Point3DTypes.h"
+
+class AllenCaloToCaloClusters final
+  : public Gaudi::Functional::MultiTransformer<std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters>(
+      const HostBuffers&)> {
+public:
+  /// Standard constructor
+  AllenCaloToCaloClusters(const std::string& name, ISvcLocator* pSvcLocator);
+
+  /// Algorithm execution
+  std::tuple<LHCb::Event::Calo::Clusters, LHCb::Event::Calo::Clusters> operator()(const HostBuffers&) const override;
+
+private:
+  Gaudi::Property<float> m_EtCalo {this, "EtCalo", 400 * Gaudi::Units::MeV, "Default ET for Calo Clusters"};
+};
+
+#endif
diff --git a/configuration/constants/Calo.json b/configuration/constants/Calo.json
new file mode 100644
index 0000000000000000000000000000000000000000..7119bb7d6c1673ae472e89be9977217e8bb946b0
--- /dev/null
+++ b/configuration/constants/Calo.json
@@ -0,0 +1,3 @@
+{
+  "configured_lines": []
+}
diff --git a/configuration/sequences/calo.py b/configuration/sequences/calo.py
new file mode 100644
index 0000000000000000000000000000000000000000..1db5eb8c1d986df120cd023b9f46b3e495952d82
--- /dev/null
+++ b/configuration/sequences/calo.py
@@ -0,0 +1,12 @@
+###############################################################################
+# (c) Copyright 2018-2020 CERN for the benefit of the LHCb Collaboration      #
+###############################################################################
+from definitions.GECSequence import GECSequence
+from definitions.CaloSequence import CaloSequence
+from definitions.algorithms import compose_sequences
+
+gec_sequence = GECSequence()
+
+calo_sequence = CaloSequence(gec_sequence['initialize_lists'])
+
+compose_sequences(gec_sequence, calo_sequence).generate()
diff --git a/configuration/sequences/definitions/CaloSequence.py b/configuration/sequences/definitions/CaloSequence.py
new file mode 100644
index 0000000000000000000000000000000000000000..0a913417627820e0e1c1d4809557b5583da8ced8
--- /dev/null
+++ b/configuration/sequences/definitions/CaloSequence.py
@@ -0,0 +1,61 @@
+###############################################################################
+# (c) Copyright 2018-2020 CERN for the benefit of the LHCb Collaboration      #
+###############################################################################
+from definitions.algorithms import *
+
+
+def CaloSequence(initialize_lists):
+    ecal_banks = data_provider_t(name="ecal_banks", bank_type="EcalPacked")
+
+    calo_count_digits = calo_count_digits_t(
+        name="calo_count_digits",
+        host_number_of_events_t=initialize_lists.host_number_of_events_t(),
+        dev_number_of_events_t=initialize_lists.dev_number_of_events_t(),
+        dev_event_list_t=initialize_lists.dev_event_list_t())
+
+    prefix_sum_ecal_num_digits = host_prefix_sum_t(
+        name="prefix_sum_ecal_num_digits",
+        dev_input_buffer_t=calo_count_digits.dev_ecal_num_digits_t())
+
+    calo_decode = calo_decode_t(
+        name="calo_decode",
+        host_number_of_events_t=initialize_lists.host_number_of_events_t(),
+        host_ecal_number_of_digits_t=prefix_sum_ecal_num_digits.
+        host_total_sum_holder_t(),
+        dev_event_list_t=initialize_lists.dev_event_list_t(),
+        dev_ecal_raw_input_t=ecal_banks.dev_raw_banks_t(),
+        dev_ecal_raw_input_offsets_t=ecal_banks.dev_raw_offsets_t(),
+        dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits.
+        dev_output_buffer_t())
+
+    calo_seed_clusters = calo_seed_clusters_t(
+        name="calo_seed_clusters",
+        host_number_of_events_t=initialize_lists.host_number_of_events_t(),
+        dev_event_list_t=initialize_lists.dev_event_list_t(),
+        dev_ecal_digits_t=calo_decode.dev_ecal_digits_t(),
+        dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits.
+        dev_output_buffer_t())
+
+    prefix_sum_ecal_num_clusters = host_prefix_sum_t(
+        name="prefix_sum_ecal_num_clusters",
+        dev_input_buffer_t=calo_seed_clusters.dev_ecal_num_clusters_t())
+
+    calo_find_clusters = calo_find_clusters_t(
+        name="calo_find_clusters",
+        host_ecal_number_of_clusters_t=prefix_sum_ecal_num_clusters.
+        host_total_sum_holder_t(),
+        dev_event_list_t=initialize_lists.dev_event_list_t(),
+        dev_ecal_digits_t=calo_decode.dev_ecal_digits_t(),
+        dev_ecal_digits_offsets_t=prefix_sum_ecal_num_digits.
+        dev_output_buffer_t(),
+        dev_ecal_seed_clusters_t=calo_seed_clusters.dev_ecal_seed_clusters_t(),
+        dev_ecal_cluster_offsets_t=prefix_sum_ecal_num_clusters.
+        dev_output_buffer_t())
+
+    calo_sequence = Sequence(ecal_banks, calo_count_digits,
+                             prefix_sum_ecal_num_digits,
+                             calo_decode,
+                             calo_seed_clusters, prefix_sum_ecal_num_clusters,
+                             calo_find_clusters)
+
+    return calo_sequence
diff --git a/device/calo/CMakeLists.txt b/device/calo/CMakeLists.txt
index a77630131e6c1008d85bda0376254c321ae581a0..32a93472a3226ec35de607871d7b9c526a3c0f65 100644
--- a/device/calo/CMakeLists.txt
+++ b/device/calo/CMakeLists.txt
@@ -9,9 +9,11 @@
 # or submit itself to any jurisdiction.                                       #
 ###############################################################################
 file(GLOB calo_decoding "decoding/src/*cu")
+file(GLOB calo_clustering "clustering/src/*cu")
 
 include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
 include_directories(${CMAKE_SOURCE_DIR}/device/event_model/calo/include)
+include_directories(clustering/include)
 include_directories(decoding/include)
 include_directories(${CMAKE_SOURCE_DIR}/main/include)
 include_directories(${CMAKE_SOURCE_DIR}/backend/include)
@@ -23,4 +25,5 @@ include_directories(${CPPGSL_INCLUDE_DIR})
 
 allen_add_device_library(Calo STATIC
   ${calo_decoding}
+  ${calo_clustering}
 )
diff --git a/device/calo/clustering/include/CaloFindClusters.cuh b/device/calo/clustering/include/CaloFindClusters.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..1d63662ede59702677f1e034c4af4603d776c95f
--- /dev/null
+++ b/device/calo/clustering/include/CaloFindClusters.cuh
@@ -0,0 +1,49 @@
+/*****************************************************************************\
+* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
+\*****************************************************************************/
+
+#pragma once
+
+#include "CaloGeometry.cuh"
+#include "CaloDigit.cuh"
+#include "CaloCluster.cuh"
+#include "DeviceAlgorithm.cuh"
+
+namespace calo_find_clusters {
+  struct Parameters {
+    HOST_INPUT(host_ecal_number_of_clusters_t, unsigned) host_ecal_number_of_clusters;
+    HOST_INPUT(host_hcal_number_of_clusters_t, unsigned) host_hcal_number_of_clusters;
+    DEVICE_INPUT(dev_event_list_t, unsigned) dev_event_list;
+    DEVICE_INPUT(dev_ecal_digits_t, CaloDigit) dev_ecal_digits;
+    DEVICE_INPUT(dev_hcal_digits_t, CaloDigit) dev_hcal_digits;
+    DEVICE_INPUT(dev_ecal_digits_offsets_t, unsigned) dev_ecal_digits_offsets;
+    DEVICE_INPUT(dev_hcal_digits_offsets_t, unsigned) dev_hcal_digits_offsets;
+    DEVICE_INPUT(dev_ecal_seed_clusters_t, CaloSeedCluster) dev_ecal_seed_clusters;
+    DEVICE_INPUT(dev_hcal_seed_clusters_t, CaloSeedCluster) dev_hcal_seed_clusters;
+    DEVICE_INPUT(dev_ecal_cluster_offsets_t, unsigned) dev_ecal_cluster_offsets;
+    DEVICE_INPUT(dev_hcal_cluster_offsets_t, unsigned) dev_hcal_cluster_offsets;
+    DEVICE_OUTPUT(dev_ecal_clusters_t, CaloCluster) dev_ecal_clusters;
+    DEVICE_OUTPUT(dev_hcal_clusters_t, CaloCluster) dev_hcal_clusters;
+    PROPERTY(block_dim_x_t, "block_dim_x", "block dimension X", unsigned) block_dim;
+  };
+
+  // Global function
+  __global__ void
+  calo_find_clusters(Parameters parameters, const char* raw_ecal_geometry, const char* raw_hcal_geometry);
+
+  // Algorithm
+  struct calo_find_clusters_t : public DeviceAlgorithm, Parameters {
+    void set_arguments_size(ArgumentReferences<Parameters>, const RuntimeOptions&, const Constants&, const HostBuffers&)
+      const;
+
+    __host__ void operator()(
+      const ArgumentReferences<Parameters>& arguments,
+      const RuntimeOptions& runtime_options,
+      const Constants& constants,
+      HostBuffers&,
+      Allen::Context const&) const;
+
+  private:
+    Property<block_dim_x_t> m_block_dim_x {this, 32};
+  };
+} // namespace calo_find_clusters
diff --git a/device/calo/clustering/include/CaloSeedClusters.cuh b/device/calo/clustering/include/CaloSeedClusters.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..cc49af82d0a2a1dd522455ce104793abe73ebc01
--- /dev/null
+++ b/device/calo/clustering/include/CaloSeedClusters.cuh
@@ -0,0 +1,58 @@
+/*****************************************************************************\
+* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
+\*****************************************************************************/
+
+#pragma once
+
+#include "CaloGeometry.cuh"
+#include "CaloDigit.cuh"
+#include "CaloCluster.cuh"
+#include "DeviceAlgorithm.cuh"
+
+namespace calo_seed_clusters {
+  struct Parameters {
+    HOST_INPUT(host_number_of_events_t, unsigned) host_number_of_events;
+    DEVICE_INPUT(dev_event_list_t, unsigned) dev_event_list;
+    DEVICE_INPUT(dev_ecal_digits_t, CaloDigit) dev_ecal_digits;
+    DEVICE_INPUT(dev_hcal_digits_t, CaloDigit) dev_hcal_digits;
+    DEVICE_INPUT(dev_ecal_digits_offsets_t, unsigned) dev_ecal_digits_offsets;
+    DEVICE_INPUT(dev_hcal_digits_offsets_t, unsigned) dev_hcal_digits_offsets;
+    DEVICE_OUTPUT(dev_ecal_num_clusters_t, unsigned) dev_ecal_num_clusters;
+    DEVICE_OUTPUT(dev_hcal_num_clusters_t, unsigned) dev_hcal_num_clusters;
+    DEVICE_OUTPUT(dev_ecal_seed_clusters_t, CaloSeedCluster) dev_ecal_seed_clusters;
+    DEVICE_OUTPUT(dev_hcal_seed_clusters_t, CaloSeedCluster) dev_hcal_seed_clusters;
+    PROPERTY(block_dim_x_t, "block_dim_x", "block dimension X", unsigned) block_dim;
+    PROPERTY(ecal_min_adc_t, "ecal_min_adc", "ECal seed cluster minimum ADC", int16_t) ecal_min_adc;
+    PROPERTY(hcal_min_adc_t, "hcal_min_adc", "HCal seed cluster minimum ADC", int16_t) hcal_min_adc;
+  };
+
+  // Global function
+  __global__ void calo_seed_clusters(
+    Parameters parameters,
+    const char* raw_ecal_geometry,
+    const char* raw_hcal_geometry,
+    const int16_t ecal_min_adc,
+    const int16_t hcal_min_adc);
+
+  // Algorithm
+  struct calo_seed_clusters_t : public DeviceAlgorithm, Parameters {
+
+    void set_arguments_size(
+      ArgumentReferences<Parameters> arguments,
+      const RuntimeOptions& runtime_options,
+      const Constants&,
+      const HostBuffers&) const;
+
+    void operator()(
+      const ArgumentReferences<Parameters>& arguments,
+      const RuntimeOptions& runtime_options,
+      const Constants& constants,
+      HostBuffers&,
+      Allen::Context const&) const;
+
+  private:
+    Property<block_dim_x_t> m_block_dim_x {this, 32};
+    Property<ecal_min_adc_t> m_ecal_min_adc {this, 20};
+    Property<hcal_min_adc_t> m_hcal_min_adc {this, 20};
+  };
+} // namespace calo_seed_clusters
diff --git a/device/calo/clustering/src/CaloFindClusters.cu b/device/calo/clustering/src/CaloFindClusters.cu
new file mode 100644
index 0000000000000000000000000000000000000000..122772215ba23d7913336638a8ccd66438b3bb7c
--- /dev/null
+++ b/device/calo/clustering/src/CaloFindClusters.cu
@@ -0,0 +1,120 @@
+/*****************************************************************************\
+* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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 <CaloCluster.cuh>
+#include <CaloFindClusters.cuh>
+
+__device__ void simple_clusters(
+  CaloDigit const* digits,
+  CaloSeedCluster const* seed_clusters,
+  CaloCluster* clusters,
+  unsigned const num_clusters,
+  const CaloGeometry& calo)
+{
+  for (unsigned c = threadIdx.x; c < num_clusters; c += blockDim.x) {
+    auto const& seed_cluster = seed_clusters[c];
+    auto& cluster = clusters[c];
+    cluster.center_id = seed_cluster.id;
+    cluster.e = calo.gain[seed_cluster.id] * (seed_cluster.adc - calo.pedestal);
+    cluster.x = seed_cluster.x;
+    cluster.y = seed_cluster.y;
+
+    uint16_t const* neighbors = &(calo.neighbors[seed_cluster.id * Calo::Constants::max_neighbours]);
+    for (uint16_t n = 0; n < Calo::Constants::max_neighbours; n++) {
+      auto const n_id = neighbors[n];
+      int16_t adc = digits[n_id].adc;
+      if (n_id != 0 && (adc != SHRT_MAX)) {
+        cluster.e += (adc - calo.pedestal) * calo.gain[n_id];
+        cluster.digits[n] = n_id;
+      }
+      else {
+        cluster.digits[n] = 0;
+      }
+    }
+
+    // Code to update position, double check
+    // for (uint16_t n = 0; n < Calo::Constants::max_neighbours; n++) {
+    //   auto const n_id = neighbors[n];
+    //   if (n_id != 0 && ((auto const adc = digits[n_id].adc) != SHRT_MAX)) {
+    //     float const adc_frac = float(adc) / float(cluster.e);
+    //     cluster.x += adc_frac * (geometry.getX(n_id) - seed_cluster.x);
+    //     cluster.y += adc_frac * (geometry.getY(n_id) - seed_cluster.y);
+    //   }
+    // }
+  }
+}
+
+__global__ void calo_find_clusters::calo_find_clusters(
+  calo_find_clusters::Parameters parameters,
+  const char* raw_ecal_geometry,
+  const char* raw_hcal_geometry)
+{
+  // Get proper geometry.
+  auto ecal_geometry = CaloGeometry(raw_ecal_geometry);
+  auto hcal_geometry = CaloGeometry(raw_hcal_geometry);
+
+  unsigned const event_number = parameters.dev_event_list[blockIdx.x];
+
+  // Build simple 3x3 clusters from seed clusters
+  // Ecal
+  unsigned const ecal_digits_offset = parameters.dev_ecal_digits_offsets[event_number];
+  unsigned const ecal_clusters_offset = parameters.dev_ecal_cluster_offsets[event_number];
+  unsigned const ecal_num_clusters = parameters.dev_ecal_cluster_offsets[event_number + 1] - ecal_clusters_offset;
+  simple_clusters(
+    parameters.dev_ecal_digits + ecal_digits_offset,
+    parameters.dev_ecal_seed_clusters + Calo::Constants::ecal_max_index / 8 * event_number,
+    parameters.dev_ecal_clusters + ecal_clusters_offset,
+    ecal_num_clusters,
+    ecal_geometry);
+
+  // Hcal
+  unsigned const hcal_digits_offset = parameters.dev_hcal_digits_offsets[event_number];
+  unsigned const hcal_clusters_offset = parameters.dev_hcal_cluster_offsets[event_number];
+  unsigned const hcal_num_clusters = parameters.dev_hcal_cluster_offsets[event_number + 1] - hcal_clusters_offset;
+  simple_clusters(
+    parameters.dev_hcal_digits + hcal_digits_offset,
+    parameters.dev_hcal_seed_clusters + Calo::Constants::hcal_max_index / 8 * event_number,
+    parameters.dev_hcal_clusters + hcal_clusters_offset,
+    hcal_num_clusters,
+    hcal_geometry);
+}
+
+void calo_find_clusters::calo_find_clusters_t::set_arguments_size(
+  ArgumentReferences<Parameters> arguments,
+  const RuntimeOptions&,
+  const Constants&,
+  const HostBuffers&) const
+{
+  set_size<dev_ecal_clusters_t>(arguments, first<host_ecal_number_of_clusters_t>(arguments));
+  set_size<dev_hcal_clusters_t>(arguments, first<host_hcal_number_of_clusters_t>(arguments));
+}
+
+__host__ void calo_find_clusters::calo_find_clusters_t::operator()(
+  const ArgumentReferences<Parameters>& arguments,
+  const RuntimeOptions& runtime_options,
+  const Constants& constants,
+  HostBuffers& host_buffers,
+  Allen::Context const& context) const
+{
+  // Find clusters.
+  initialize<dev_ecal_clusters_t>(arguments, SHRT_MAX, context);
+  initialize<dev_hcal_clusters_t>(arguments, SHRT_MAX, context);
+
+  global_function(calo_find_clusters)(
+    dim3(size<dev_event_list_t>(arguments)), dim3(property<block_dim_x_t>().get()), context)(
+    arguments, constants.dev_ecal_geometry, constants.dev_hcal_geometry);
+
+  if (runtime_options.do_check) {
+    safe_assign_to_host_buffer<dev_ecal_cluster_offsets_t>(host_buffers.host_ecal_cluster_offsets, arguments, context);
+    safe_assign_to_host_buffer<dev_hcal_cluster_offsets_t>(host_buffers.host_hcal_cluster_offsets, arguments, context);
+    safe_assign_to_host_buffer<dev_ecal_clusters_t>(host_buffers.host_ecal_clusters, arguments, context);
+    safe_assign_to_host_buffer<dev_hcal_clusters_t>(host_buffers.host_hcal_clusters, arguments, context);
+  }
+}
diff --git a/device/calo/clustering/src/CaloSeedClusters.cu b/device/calo/clustering/src/CaloSeedClusters.cu
new file mode 100644
index 0000000000000000000000000000000000000000..d145a3df8f6c6124b51a2c023c2691b46c711713
--- /dev/null
+++ b/device/calo/clustering/src/CaloSeedClusters.cu
@@ -0,0 +1,107 @@
+/*****************************************************************************\
+* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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 <CaloSeedClusters.cuh>
+
+__device__ void seed_clusters(
+  CaloDigit const* digits,
+  unsigned const num_digits,
+  CaloSeedCluster* clusters,
+  unsigned* num_clusters,
+  const CaloGeometry& geometry,
+  const int16_t min_adc)
+{
+  // Loop over all CellIDs.
+  for (unsigned i = threadIdx.x; i < num_digits; i += blockDim.x) {
+    int16_t adc = digits[i].adc;
+    if (adc == SHRT_MAX || adc < min_adc) {
+      continue;
+    }
+    uint16_t* neighbors = &(geometry.neighbors[i * Calo::Constants::max_neighbours]);
+    bool is_max = true;
+    for (unsigned n = 0; n < Calo::Constants::max_neighbours; n++) {
+      auto const neighbor_adc = digits[neighbors[n]].adc;
+      is_max = is_max && (neighbor_adc == SHRT_MAX || adc > neighbor_adc);
+    }
+    if (is_max) {
+      auto const id = atomicAdd(num_clusters, 1);
+      clusters[id] = CaloSeedCluster(i, digits[i].adc, geometry.getX(i), geometry.getY(i));
+    }
+  }
+}
+
+__global__ void calo_seed_clusters::calo_seed_clusters(
+  calo_seed_clusters::Parameters parameters,
+  const char* raw_ecal_geometry,
+  const char* raw_hcal_geometry,
+  const int16_t ecal_min_adc,
+  const int16_t hcal_min_adc)
+{
+  unsigned const event_number = parameters.dev_event_list[blockIdx.x];
+
+  // Get geometry.
+  auto ecal_geometry = CaloGeometry(raw_ecal_geometry);
+  auto hcal_geometry = CaloGeometry(raw_hcal_geometry);
+
+  // ECal
+  auto const ecal_digits_offset = parameters.dev_ecal_digits_offsets[event_number];
+  seed_clusters(
+    &parameters.dev_ecal_digits[ecal_digits_offset],
+    parameters.dev_ecal_digits_offsets[event_number + 1] - ecal_digits_offset,
+    &parameters.dev_ecal_seed_clusters[Calo::Constants::ecal_max_index / 8 * event_number],
+    &parameters.dev_ecal_num_clusters[event_number],
+    ecal_geometry,
+    ecal_min_adc);
+
+  // HCal
+  auto const hcal_digits_offset = parameters.dev_hcal_digits_offsets[event_number];
+  seed_clusters(
+    &parameters.dev_hcal_digits[hcal_digits_offset],
+    parameters.dev_hcal_digits_offsets[event_number + 1] - hcal_digits_offset,
+    &parameters.dev_hcal_seed_clusters[Calo::Constants::hcal_max_index / 8 * event_number],
+    &parameters.dev_hcal_num_clusters[event_number],
+    hcal_geometry,
+    hcal_min_adc);
+}
+
+void calo_seed_clusters::calo_seed_clusters_t::set_arguments_size(
+  ArgumentReferences<Parameters> arguments,
+  const RuntimeOptions&,
+  const Constants&,
+  const HostBuffers&) const
+{
+  auto const n_events = first<host_number_of_events_t>(arguments);
+  set_size<dev_ecal_num_clusters_t>(arguments, n_events);
+  set_size<dev_hcal_num_clusters_t>(arguments, n_events);
+
+  // TODO: get this from the geometry too
+  set_size<dev_ecal_seed_clusters_t>(arguments, Calo::Constants::ecal_max_index / 8 * n_events);
+  set_size<dev_hcal_seed_clusters_t>(arguments, Calo::Constants::hcal_max_index / 8 * n_events);
+}
+
+void calo_seed_clusters::calo_seed_clusters_t::operator()(
+  const ArgumentReferences<Parameters>& arguments,
+  const RuntimeOptions&,
+  const Constants& constants,
+  HostBuffers&,
+  Allen::Context const& context) const
+{
+  initialize<dev_ecal_num_clusters_t>(arguments, 0, context);
+  initialize<dev_hcal_num_clusters_t>(arguments, 0, context);
+
+  // Find local maxima.
+  global_function(calo_seed_clusters)(
+    dim3(size<dev_event_list_t>(arguments)), dim3(property<block_dim_x_t>().get()), context)(
+    arguments,
+    constants.dev_ecal_geometry,
+    constants.dev_hcal_geometry,
+    property<ecal_min_adc_t>().get(),
+    property<hcal_min_adc_t>().get());
+}
diff --git a/device/event_model/calo/include/CaloCluster.cuh b/device/event_model/calo/include/CaloCluster.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..d0bc0c09c644a547da0ab18eae7683b599621464
--- /dev/null
+++ b/device/event_model/calo/include/CaloCluster.cuh
@@ -0,0 +1,49 @@
+/*****************************************************************************\
+* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
+*                                                                             *
+* This software is distributed under the terms of the Apache License          *
+* version 2 (Apache-2.0), copied verbatim in the file "LICENSE".              *
+*                                                                             *
+* 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.                                       *
+\*****************************************************************************/
+#pragma once
+
+#include <BackendCommon.h>
+#include "CaloConstants.cuh"
+
+struct CaloDigitClusters {
+  uint16_t clustered_at_iteration = 0;
+  uint16_t clusters[Calo::Constants::digit_max_clusters];
+
+  __device__ __host__ CaloDigitClusters()
+  {
+    for (int i = 0; i < Calo::Constants::digit_max_clusters; i++) {
+      clusters[i] = 0;
+    }
+  }
+};
+
+struct CaloCluster {
+  float e = 0.f;
+  float x = 0.f;
+  float y = 0.f;
+  uint16_t center_id = 0;
+  uint16_t digits[Calo::Constants::max_neighbours] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+  __device__ __host__ CaloCluster(uint16_t index, float energy, float rX, float rY) :
+    e {energy}, x {rX}, y {rY}, center_id {index}
+  {}
+};
+
+struct CaloSeedCluster {
+  uint16_t id = 0;
+  int16_t adc = 0; // Is a double in the original algorithm, but is still integer here?
+  float x = 0.f;
+  float y = 0.f;
+
+  __device__ __host__ CaloSeedCluster(uint16_t cellid, int16_t a, float rX, float rY) :
+    id {cellid}, adc {a}, x {rX}, y {rY}
+  {}
+};
diff --git a/device/event_model/calo/include/CaloConstants.cuh b/device/event_model/calo/include/CaloConstants.cuh
index f43e134c7654abb762721cbc96e0aca444de696a..b16d2bbcbcaba7c6436e3277f8f630c950eb6a5c 100644
--- a/device/event_model/calo/include/CaloConstants.cuh
+++ b/device/event_model/calo/include/CaloConstants.cuh
@@ -13,5 +13,16 @@
 namespace Calo {
   namespace Constants {
     constexpr uint16_t card_channels = 32;
+    constexpr uint16_t ecal_max_index = 6016;
+    constexpr uint16_t hcal_max_index = 1488;
+
+    // Max distance based on CellIDs is 64 steps away, so the iteration in which a cell is clustered can never be more
+    // than 64.
+    constexpr uint16_t unclustered = 65;
+    constexpr uint16_t ignore = unclustered + 1;
+
+    constexpr uint16_t digit_max_clusters = 15;
+    constexpr uint16_t max_neighbours = 9;
+
   } // namespace Constants
 } // namespace Calo
diff --git a/main/src/Allen.cpp b/main/src/Allen.cpp
index c7d589db99978b5c623913cea5efcda05329cde2..ad5290ce01b6b56033ce889a9b090622737aecd0 100644
--- a/main/src/Allen.cpp
+++ b/main/src/Allen.cpp
@@ -388,8 +388,13 @@ int allen(
   else {
     mep_layout = false;
     // The binary input provider expects the folders for the bank types as connections
-    std::vector<std::string> connections = {
-      folder_name_velopix_raw, folder_name_UT_raw, folder_name_SciFi_raw, folder_name_Muon_raw, folder_name_ODIN_raw};
+    std::vector<std::string> connections = {folder_name_velopix_raw,
+                                            folder_name_UT_raw,
+                                            folder_name_SciFi_raw,
+                                            folder_name_Muon_raw,
+                                            folder_name_ecal_raw,
+                                            folder_name_hcal_raw,
+                                            folder_name_ODIN_raw};
     input_provider =
       std::make_unique<BinaryProvider<BankTypes::VP, BankTypes::UT, BankTypes::FT, BankTypes::MUON, BankTypes::ODIN>>(
         number_of_slices, events_per_slice, n_events, std::move(connections), n_io_reps, file_list);
diff --git a/stream/sequence/include/HostBuffers.cuh b/stream/sequence/include/HostBuffers.cuh
index 7cb456619e49787111e157dd4361989f274c3d50..7d06b34757002383fa870857881bfae71b0c0f9d 100644
--- a/stream/sequence/include/HostBuffers.cuh
+++ b/stream/sequence/include/HostBuffers.cuh
@@ -10,6 +10,7 @@
 #include <gsl/gsl>
 
 #include <CaloDigit.cuh>
+#include <CaloCluster.cuh>
 
 #include <BackendCommon.h>
 
@@ -109,6 +110,10 @@ struct HostBuffers {
   bool* host_match_upstream_muon;
 
   // Calo
+  gsl::span<unsigned> host_ecal_cluster_offsets = {};
+  gsl::span<unsigned> host_hcal_cluster_offsets = {};
+  gsl::span<CaloCluster> host_ecal_clusters = {};
+  gsl::span<CaloCluster> host_hcal_clusters = {};
   gsl::span<unsigned> host_ecal_digits_offsets = {};
   gsl::span<unsigned> host_hcal_digits_offsets = {};
   gsl::span<CaloDigit> host_ecal_digits = {};