diff --git a/Analysis/CMakeLists.txt b/Analysis/CMakeLists.txt
index 9009dbfc5d635a07ba3d53e03aa28b8ca2b2af06..4973572604d41cdae50f054810e456730ab81a0c 100644
--- a/Analysis/CMakeLists.txt
+++ b/Analysis/CMakeLists.txt
@@ -3,3 +3,4 @@ add_subdirectory(TotalInvMass)
 add_subdirectory(TrackInspect)
 add_subdirectory(DumpEvent)
 add_subdirectory(ReadDigi)
+add_subdirectory(JetClustering)
diff --git a/Analysis/JetClustering/CMakeLists.txt b/Analysis/JetClustering/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d4df15822747f001b4318d05384769a7119f1d0b
--- /dev/null
+++ b/Analysis/JetClustering/CMakeLists.txt
@@ -0,0 +1,21 @@
+# Set the CMake module path if necessary
+
+# Modules
+gaudi_add_module(JetClustering
+                 SOURCES src/JetClustering.cpp
+                 LINK GearSvc
+                      Gaudi::GaudiKernel
+                      k4FWCore::k4FWCore
+                      DetInterface
+                      DataHelperLib
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${ROOT_LIBRARIES}
+                      ${FASTJET_LIBRARY}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)
+install(TARGETS JetClustering
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)
diff --git a/Analysis/JetClustering/src/JetClustering.cpp b/Analysis/JetClustering/src/JetClustering.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3c2c8133040402a5eeaea28589104c8843b771a
--- /dev/null
+++ b/Analysis/JetClustering/src/JetClustering.cpp
@@ -0,0 +1,190 @@
+#include "JetClustering.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+#include "DetInterface/IGeomSvc.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/PseudoJet.hh"
+//time measurement
+#include <chrono>
+
+using namespace fastjet;
+using namespace std;
+
+DECLARE_COMPONENT( JetClustering )
+
+//------------------------------------------------------------------------------
+JetClustering::JetClustering( const string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+  declareProperty("InputPFOs", m_PFOColHdl, "Handle of the Input PFO collection");
+  declareProperty("Algorithm", m_algo, "Jet clustering algorithm");
+  declareProperty("nJets", m_nJets, "Number of jets to be clustered");
+  declareProperty("R", m_R, "Jet clustering radius");
+  declareProperty("OutputFile", m_outputFile, "Output file name");
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::initialize(){
+
+  _nEvt = 0;
+
+  // Create a new TTree
+  _file = TFile::Open(m_outputFile.value().c_str(), "RECREATE");
+  _tree = new TTree("jets", "jets");
+  _tree->Branch("jet1_px", &jet1_px, "jet1_px/D");
+  _tree->Branch("jet1_py", &jet1_py, "jet1_py/D");
+  _tree->Branch("jet1_pz", &jet1_pz, "jet1_pz/D");
+  _tree->Branch("jet1_E", &jet1_E, "jet1_E/D");
+  _tree->Branch("jet1_costheta", &jet1_costheta, "jet1_costheta/D");
+  _tree->Branch("jet1_phi", &jet1_phi, "jet1_phi/D");
+  _tree->Branch("jet1_pt", &jet1_pt, "jet1_pt/D");
+  _tree->Branch("jet1_nconstituents", &jet1_nconstituents, "jet1_nconstituents/I");
+  _tree->Branch("jet2_px", &jet2_px, "jet2_px/D");
+  _tree->Branch("jet2_py", &jet2_py, "jet2_py/D");
+  _tree->Branch("jet2_pz", &jet2_pz, "jet2_pz/D");
+  _tree->Branch("jet2_E", &jet2_E, "jet2_E/D");
+  _tree->Branch("jet2_costheta", &jet2_costheta, "jet2_costheta/D");
+  _tree->Branch("jet2_phi", &jet2_phi, "jet2_phi/D");
+  _tree->Branch("jet2_pt", &jet2_pt, "jet2_pt/D");
+  _tree->Branch("jet2_nconstituents", &jet2_nconstituents, "jet2_nconstituents/I");
+  _tree->Branch("constituents_E1tot", &constituents_E1tot, "constituents_E1tot/D");
+  _tree->Branch("constituents_E2tot", &constituents_E2tot, "constituents_E2tot/D");
+  _tree->Branch("mass", &mass, "mass/D");
+  _tree->Branch("ymerge", ymerge, "ymerge[6]/D");
+  //_tree->Branch("px", &px);
+  //_tree->Branch("py", &py);
+  //_tree->Branch("pz", &pz);
+  //_tree->Branch("E", &E);
+  
+  return StatusCode::SUCCESS;
+
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::execute(){
+  
+  const edm4hep::ReconstructedParticleCollection* PFO = nullptr;
+  try {
+    PFO = m_PFOColHdl.get();
+  }
+  catch ( GaudiException &e ){
+    info() << "Collection " << m_PFOColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  // Check if the collection is empty
+  if ( PFO->size() == 0 ) {
+    info() << "Collection " << m_PFOColHdl.fullKey() << " is empty in event " << _nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  //initial indces
+  CleanVars();
+
+  //resize particles. there are many particles are empty, maybe upstream bugs
+  vector<PseudoJet> input_particles;
+  for(const auto& pfo : *PFO){
+    if ( isnan( pfo.getMomentum()[0]) || isnan(pfo.getMomentum()[1]) || isnan(pfo.getMomentum()[2]) || pfo.getEnergy() == 0 ) {
+      debug() << "Particle " << _particle_index << " px: " << pfo.getMomentum()[0] << " py: " << pfo.getMomentum()[1] << " pz: " << pfo.getMomentum()[2] << " E: " << pfo.getEnergy()  << endmsg;
+      _particle_index++;
+      continue;
+    }
+    else{
+      input_particles.push_back( PseudoJet( pfo.getMomentum()[0], pfo.getMomentum()[1], pfo.getMomentum()[2], pfo.getEnergy() ) );
+      _particle_index++;
+      debug() << "Particle " << _particle_index << " px: "<<pfo.getMomentum()[0]<<" py: "<<pfo.getMomentum()[1]<<" pz: "<<pfo.getMomentum()[2]<<" E: "<<pfo.getEnergy()<<endmsg;
+    }
+  }
+
+  // create a jet definition
+  int nJets = m_nJets;
+  double R = m_R;
+
+  JetDefinition jet_def(ee_kt_algorithm);
+  //JetDefinition jet_def(antikt_algorithm, R);
+
+  // run the clustering, extract the jets
+  //std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
+
+  ClusterSequence clust_seq(input_particles, jet_def);
+  
+  vector<PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(nJets));
+  
+  //std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
+  
+  //_time = std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1000000.0;
+  
+  //info() << "FastJet clustering done in seconds: " << _time << endmsg;
+  
+  //string jet_def_str = jet_def.description();
+  //debug() << "Clustering with " << jet_def_str << " for " << nJets << " jets with R = " << R << endmsg;
+  
+  // print the jets
+  debug() << "                     px           py           pz           E           costheta  phi " << endmsg;
+  for (unsigned i = 0; i < jets.size(); i++) {
+    
+    debug() << "   jet " << i << "             " << jets[i].px() << " " << jets[i].py() << " " << jets[i].pz() << " " << jets[i].E() << " " << jets[i].cos_theta() << " " << jets[i].phi() << endmsg;
+    
+    vector<PseudoJet> constituents = jets[i].constituents();
+    
+    for (unsigned j = 0; j < constituents.size(); j++) {
+      
+      if ( i == 0 )constituents_E1tot += constituents[j].E();
+      
+      if ( i == 1 )constituents_E2tot += constituents[j].E();
+      
+      debug() << "    constituent  " << j << " " << constituents[j].px() << " " << constituents[j].py() << " " << constituents[j].pz() << " " << constituents[j].E() << " " << constituents[j].cos_theta() << " " << constituents[j].phi() << endmsg;
+    
+    }
+  
+  }
+  
+  debug() << "Total energy of constituents: " << constituents_E1tot << " " << constituents_E2tot << endmsg;  
+
+  double _ymin[20];
+  for(int i=1; i<6;i++){
+    _ymin[i-1] = clust_seq.exclusive_ymerge (i);
+    debug() << " -log10(y" << i << i+1 << ") = " << -log10(_ymin[i-1]) << endmsg;
+  }
+  mass = (jets[0] + jets[1]).m();
+
+  info() << "Total mass of all particles: " << " mass of jet1+jet2: " << mass << endmsg;
+
+  // fill the tree
+  jet1_px = jets[0].px();
+  jet1_py = jets[0].py();
+  jet1_pz = jets[0].pz();
+  jet1_E = jets[0].E();
+  jet1_costheta = jets[0].cos_theta();
+  jet1_phi = jets[0].phi();
+  jet1_pt = jets[0].pt();
+  jet1_nconstituents = jets[0].constituents().size();
+  jet2_px = jets[1].px();
+  jet2_py = jets[1].py();
+  jet2_pz = jets[1].pz();
+  jet2_E = jets[1].E();
+  jet2_costheta = jets[1].cos_theta();
+  jet2_phi = jets[1].phi();
+  jet2_pt = jets[1].pt();
+  jet2_nconstituents = jets[1].constituents().size();
+  for(int i=0; i<6; i++){
+    ymerge[i] = _ymin[i];
+  }
+  _tree->Fill();
+
+  _nEvt++;
+  return StatusCode::SUCCESS;
+
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::finalize(){
+  debug() << "Finalizing..." << endmsg;
+  _file->Write();
+  return StatusCode::SUCCESS;
+}
diff --git a/Analysis/JetClustering/src/JetClustering.h b/Analysis/JetClustering/src/JetClustering.h
new file mode 100644
index 0000000000000000000000000000000000000000..5d0672909eee7cac495a6318bdce25182638ebe3
--- /dev/null
+++ b/Analysis/JetClustering/src/JetClustering.h
@@ -0,0 +1,75 @@
+#ifndef JetClustering_h
+#define JetClustering_h 1
+
+#include "UTIL/ILDConf.h"
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/Algorithm.h"
+#include <random>
+#include "GaudiKernel/NTuple.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include "edm4hep/ReconstructedParticleData.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "edm4hep/ReconstructedParticle.h"
+
+class JetClustering : public Algorithm {
+ public:
+  // Constructor of this form must be provided
+  JetClustering( const std::string& name, ISvcLocator* pSvcLocator );
+
+  // Three mandatory member functions of any algorithm
+  StatusCode initialize() override;
+  StatusCode execute() override;
+  StatusCode finalize() override;
+
+ private:
+  DataHandle<edm4hep::ReconstructedParticleCollection> m_PFOColHdl{"PandoraPFOs", Gaudi::DataHandle::Reader, this};
+  Gaudi::Property<std::string> m_algo{this, "Algorithm", "ee_kt_algorithm"};
+  Gaudi::Property<int> m_nJets{this, "nJets", 2};
+  Gaudi::Property<double> m_R{this, "R", 0.6};
+  Gaudi::Property<std::string> m_outputFile{this, "OutputFile", "JetClustering.root"};
+
+  int _nEvt;
+  int _particle_index;
+  int _jet_index;
+  int _nparticles;
+  double _time;
+
+  TFile* _file;
+  TTree* _tree;
+  double jet1_px, jet1_py, jet1_pz, jet1_E;
+  double jet2_px, jet2_py, jet2_pz, jet2_E;
+  double jet1_costheta, jet1_phi;
+  double jet2_costheta, jet2_phi;
+  double jet1_pt, jet2_pt;
+  int jet1_nconstituents, jet2_nconstituents;
+  double constituents_E1tot;
+  double constituents_E2tot;
+  double ymerge[6];
+  double mass;
+
+  void CleanVars(){
+
+    _particle_index = 0;
+    _jet_index = 0;
+    _nparticles = 0;
+    _time = 0;
+    
+    jet1_px = jet1_py = jet1_pz = jet1_E = 0;
+    jet2_px = jet2_py = jet2_pz = jet2_E = 0;
+    jet1_costheta = jet1_phi = 0;
+    jet2_costheta = jet2_phi = 0;
+    jet1_pt = jet2_pt = 0;
+    jet1_nconstituents = jet2_nconstituents = 0;
+    constituents_E1tot = constituents_E2tot = 0;
+    for(int i=0; i<6; i++){
+      ymerge[i] = 0;
+    }
+    mass = 0;
+
+  }
+
+};
+#endif
+
diff --git a/Examples/options/jcls.py b/Examples/options/jcls.py
new file mode 100644
index 0000000000000000000000000000000000000000..c20c7005cb7498006fd5eb9ad9ba1ce6b542c29d
--- /dev/null
+++ b/Examples/options/jcls.py
@@ -0,0 +1,34 @@
+import os, sys
+from Gaudi.Configuration import *
+
+########### k4DataSvc ####################
+from Configurables import k4DataSvc
+podioevent = k4DataSvc("EventDataSvc", input="Rec_TDR_o1_v01_E240_nnHgg.root")
+##########################################
+
+########## CEPCSWData ################# 
+cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
+#######################################
+
+########## Podio Input ###################
+from Configurables import PodioInput
+inp = PodioInput("InputReader")
+inp.collections = [ "PandoraPFOs" ]
+##########################################
+
+from Configurables import JetClustering
+jetclustering = JetClustering("JetClustering")
+jetclustering.nJets = 2
+jetclustering.R = 0.6
+jetclustering.OutputFile = "Rec_JetClustering_TDR_o1_v01_E240_nnHgg.root"
+
+########################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( 
+    TopAlg=[inp, jetclustering ],
+    EvtSel="NONE",
+    EvtMax=100,
+    ExtSvc=[podioevent],
+    #OutputLevel=DEBUG
+)
\ No newline at end of file
diff --git a/cmake/CEPCSWDependencies.cmake b/cmake/CEPCSWDependencies.cmake
index ed2b35ccc65f61c484aa9bcf59ac38aac8a89310..b8d28b7d7a782a6b1de884a13cb405b6af55e86d 100644
--- a/cmake/CEPCSWDependencies.cmake
+++ b/cmake/CEPCSWDependencies.cmake
@@ -62,3 +62,5 @@ else()
   message("Try to use an internal installation of EDM4CEPC")
   include("${CMAKE_CURRENT_LIST_DIR}/internal_edm4cepc.cmake")
 endif()
+
+find_package(FastJet)
diff --git a/cmake/FindFastJet.cmake b/cmake/FindFastJet.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..16ce192365e54f6f30e56aedd90b862928c24744
--- /dev/null
+++ b/cmake/FindFastJet.cmake
@@ -0,0 +1,29 @@
+# Source code from:
+#   https://github.com/HSF/cmaketools/blob/master/modules/FindFastJet.cmake
+#
+# - Locate FastJet library
+# Defines:
+#
+#  FASTJET_FOUND
+#  FASTJET_INCLUDE_DIR
+#  FASTJET_INCLUDE_DIRS (not cached)
+#  FASTJET_LIBRARY
+#  FASTJET_LIBRARIES (not cached)
+#  FASTJET_LIBRARY_DIRS (not cached)
+
+find_path(FASTJET_INCLUDE_DIR fastjet/version.hh
+          HINTS $ENV{FASTJET_ROOT_DIR}/include ${FASTJET_ROOT_DIR}/include)
+
+find_library(FASTJET_LIBRARY NAMES fastjet
+             HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib)
+                   
+# handle the QUIETLY and REQUIRED arguments and set FASTJET_FOUND to TRUE if
+# all listed variables are TRUE
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(FastJet DEFAULT_MSG FASTJET_INCLUDE_DIR FASTJET_LIBRARY)
+
+mark_as_advanced(FASTJET_FOUND FASTJET_INCLUDE_DIR FASTJET_LIBRARY)
+
+set(FASTJET_INCLUDE_DIRS ${FASTJET_INCLUDE_DIR})
+set(FASTJET_LIBRARIES ${FASTJET_LIBRARY})
+get_filename_component(FASTJET_LIBRARY_DIRS ${FASTJET_LIBRARY} PATH)
\ No newline at end of file