diff --git a/.gitignore b/.gitignore
index 58e07a1f2..23b5389d2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@ build
 spack*
 ./Generator/output/
 ./Generator/options/
+InstallArea/
diff --git a/Digitisers/CMakeLists.txt b/Digitisers/CMakeLists.txt
index 6f6f2e404..44b4b03ba 100644
--- a/Digitisers/CMakeLists.txt
+++ b/Digitisers/CMakeLists.txt
@@ -2,3 +2,4 @@ add_subdirectory(DCHDigi)
 add_subdirectory(G2CDArbor)
 add_subdirectory(SimHitMerge)
 add_subdirectory(SimpleDigi)
+add_subdirectory(CRDEcalDigi)
diff --git a/Digitisers/CRDEcalDigi/CMakeLists.txt b/Digitisers/CRDEcalDigi/CMakeLists.txt
new file mode 100644
index 000000000..4a44558e2
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/CMakeLists.txt
@@ -0,0 +1,22 @@
+# Modules
+gaudi_add_module(CRDEcalDigi
+                 SOURCES src/CRDEcalDigiAlg.cpp
+                 SOURCES src/CRDHcalDigiAlg.cpp
+                 LINK k4FWCore::k4FWCore
+                      GearSvc
+                      DetInterface
+                      Gaudi::GaudiKernel
+                      Gaudi::GaudiAlgLib 
+                      ${CLHEP_LIBRARIES}
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${LCIO_LIBRARIES}
+                      ${ROOT_LIBRARIES}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)
+install(TARGETS CRDEcalDigi
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)
+
diff --git a/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp
new file mode 100644
index 000000000..f52c1153b
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp
@@ -0,0 +1,505 @@
+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+// Unit in code: mm, ns. 
+// NOTE: This digitialization highly matches detector geometry CRDEcalBarrel_v01.
+// TODO: read geometry info automatically.  
+#include "CRDEcalDigiAlg.h"
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+#include "edm4hep/Cluster.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+
+#include "TVector3.h"
+#include <math.h>
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <map>
+
+// #include <fstream>
+// #include <ctime>
+
+#define C 299.79  // unit: mm/ns
+#define PI 3.141592653
+using namespace std;
+using namespace dd4hep;
+
+DECLARE_COMPONENT( CRDEcalDigiAlg )
+
+CRDEcalDigiAlg::CRDEcalDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+  : GaudiAlgorithm(name, svcLoc),
+    _nEvt(0)
+{
+  
+	// Input collections
+	declareProperty("SimCaloHitCollection", r_SimCaloCol, "Handle of the Input SimCaloHit collection");
+  
+	// Output collections
+	declareProperty("CaloHitCollection", w_DigiCaloCol, "Handle of Digi CaloHit collection");
+	declareProperty("CaloAssociationCollection", w_CaloAssociationCol, "Handle of CaloAssociation collection");
+	declareProperty("CaloMCPAssociationCollection", w_MCPCaloAssociationCol, "Handle of CaloAssociation collection");
+   
+}
+
+StatusCode CRDEcalDigiAlg::initialize()
+{
+
+  if(_writeNtuple){
+    std::string s_outfile = _filename;
+    m_wfile = new TFile(s_outfile.c_str(), "recreate");
+    t_SimCont = new TTree("SimStep", "SimStep");
+    t_SimBar = new TTree("SimBarHit", "SimBarHit");
+    t_SimCont->Branch("step_x", &m_step_x);
+    t_SimCont->Branch("step_y", &m_step_y);
+    t_SimCont->Branch("step_z", &m_step_z);
+    t_SimCont->Branch("step_t", &m_step_t);			// yyy: time of each step
+    t_SimCont->Branch("stepBar_x", &m_stepBar_x);
+    t_SimCont->Branch("stepBar_y", &m_stepBar_y);
+    t_SimCont->Branch("stepBar_z", &m_stepBar_z);
+    t_SimCont->Branch("step_E", &m_step_E);
+    t_SimCont->Branch("step_T1", &m_step_T1);
+    t_SimCont->Branch("step_T2", &m_step_T2);
+    t_SimBar->Branch("totE", &totE);
+    t_SimBar->Branch("simBar_x", &m_simBar_x);
+    t_SimBar->Branch("simBar_y", &m_simBar_y);
+    t_SimBar->Branch("simBar_z", &m_simBar_z);
+    t_SimBar->Branch("simBar_T1", &m_simBar_T1);
+    t_SimBar->Branch("simBar_T2", &m_simBar_T2);
+    t_SimBar->Branch("simBar_Q1", &m_simBar_Q1);
+    t_SimBar->Branch("simBar_Q2", &m_simBar_Q2);
+    t_SimBar->Branch("simBar_module", &m_simBar_module);
+    t_SimBar->Branch("simBar_stave", &m_simBar_stave);
+    t_SimBar->Branch("simBar_dlayer", &m_simBar_dlayer);
+    t_SimBar->Branch("simBar_part", &m_simBar_part);
+    t_SimBar->Branch("simBar_slayer", &m_simBar_slayer);
+    t_SimBar->Branch("simBar_cellID", &m_simBar_cellID);
+  }
+
+	std::cout<<"CRDEcalDigiAlg::m_scale="<<m_scale<<std::endl;
+	m_geosvc = service<IGeomSvc>("GeomSvc");
+	if ( !m_geosvc )  throw "CRDEcalDigiAlg :Failed to find GeomSvc ...";
+	dd4hep::Detector* m_dd4hep = m_geosvc->lcdd();
+	if ( !m_dd4hep )  throw "CRDEcalDigiAlg :Failed to get dd4hep::Detector ...";
+	m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep);
+   m_decoder = m_geosvc->getDecoder(_readoutName);
+	if (!m_decoder) {
+		error() << "Failed to get the decoder. " << endmsg;
+		return StatusCode::FAILURE;
+	}
+
+	//m_edmsvc = service<ICRDEcalSvc>("CRDEcalSvc");
+	//if ( !m_edmsvc )  throw "CRDEcalDigiAlg :Failed to find CRDEcalSvc ...";
+
+	rndm.SetSeed(_seed);
+	std::cout<<"CRDEcalDigiAlg::initialize"<<std::endl;
+	return GaudiAlgorithm::initialize();
+}
+
+StatusCode CRDEcalDigiAlg::execute()
+{
+// clock_t yyy_start, yyy_enddigi;
+// yyy_start = clock(); // 记录开始时间
+
+	if(_nEvt==0) std::cout<<"CRDEcalDigiAlg::execute Start"<<std::endl;
+	std::cout<<"Processing event: "<<_nEvt<<std::endl;
+   if(_nEvt<_Nskip){ _nEvt++; return StatusCode::SUCCESS; }
+
+	Clear();
+
+ 	const edm4hep::SimCalorimeterHitCollection* SimHitCol =  r_SimCaloCol.get();
+
+	edm4hep::CalorimeterHitCollection* caloVec = w_DigiCaloCol.createAndPut();
+	edm4hep::MCRecoCaloAssociationCollection* caloAssoVec = w_CaloAssociationCol.createAndPut();
+  edm4hep::MCRecoCaloParticleAssociationCollection* caloMCPAssoVec = w_MCPCaloAssociationCol.createAndPut();
+ 	std::vector<edm4hep::SimCalorimeterHit> m_simhitCol; m_simhitCol.clear();
+  std::vector<CaloBar> m_barCol; m_barCol.clear(); 
+
+  if(SimHitCol == 0) 
+  {
+     std::cout<<"not found SimCalorimeterHitCollection"<< std::endl;
+     return StatusCode::SUCCESS;
+  }
+  if(_Debug>=1) std::cout<<"digi, input sim hit size="<< SimHitCol->size() <<std::endl;
+
+	totE=0;
+	double totE_Digi=0;
+
+	//Merge input simhit(steps) to real simhit(bar).
+	MergeHits(*SimHitCol, m_simhitCol);
+	if(_Debug>=1) std::cout<<"Finish Hit Merge, with Nhit: "<<m_simhitCol.size()<<std::endl;
+
+	//Loop in SimHit, digitalize SimHit to DigiBar
+	for(int i=0;i<m_simhitCol.size();i++){
+
+		auto SimHit = m_simhitCol.at(i);
+		if(SimHit.getEnergy()<_Eth) continue;
+
+		unsigned long long id = SimHit.getCellID();
+		CaloBar hitbar;
+		hitbar.setcellID( id);
+		hitbar.setcellID(	m_decoder->get(id, "system"), 
+											m_decoder->get(id, "module"), 
+											m_decoder->get(id, "stave"), 
+											m_decoder->get(id, "dlayer"), 
+											m_decoder->get(id, "part"), 
+											m_decoder->get(id, "slayer"),
+											m_decoder->get(id, "bar"));
+
+		double Lbar = GetBarLength(hitbar);  //NOTE: Is fixed with geometry CRDEcalBarrel_v01. 
+		dd4hep::Position hitpos = m_cellIDConverter->position(id);
+    TVector3 barpos(10*hitpos.x(), 10*hitpos.y(), 10*hitpos.z()); //cm to mm.
+		hitbar.setPosition(barpos);
+		//hitbar.T1 = 99999; hitbar.T2 = 99999;
+		//if(_Debug>=2) std::cout<<"SimHit contribution size: "<<SimHit.contributions_size()<<std::endl;
+
+    MCParticleToEnergyWeightMap MCPEnMap; MCPEnMap.clear();
+		std::vector<HitStep> DigiLvec; DigiLvec.clear();
+		std::vector<HitStep> DigiRvec; DigiRvec.clear();
+		double totQ1 = 0;
+		double totQ2 = 0;
+
+		//Loop in all SimHitContribution(G4Step). 
+		for(int iCont=0; iCont < SimHit.contributions_size(); ++iCont){
+			auto conb = SimHit.getContributions(iCont);
+			if( !conb.isAvailable() ) { std::cout<<"CRDEcalDigiAlg  Can not get SimHitContribution: "<<iCont<<std::endl; continue;}
+
+			double en = conb.getEnergy();
+			if(en == 0) continue;
+
+      auto mcp = conb.getParticle();
+      MCPEnMap[mcp] += en;
+			TVector3 steppos(conb.getStepPosition().x, conb.getStepPosition().y, conb.getStepPosition().z);
+			TVector3 rpos = steppos-hitbar.getPosition();
+			float step_time = conb.getTime();		// yyy: step time
+
+			m_step_x.push_back(steppos.x());
+			m_step_y.push_back(steppos.y());
+			m_step_z.push_back(steppos.z());
+			m_step_t.push_back(step_time);			// yyy: push back step time
+			m_step_E.push_back(en);
+			m_stepBar_x.push_back(hitbar.getPosition().x());
+			m_stepBar_y.push_back(hitbar.getPosition().y());
+			m_stepBar_z.push_back(hitbar.getPosition().z());
+
+			if(_Debug>=3){
+				cout<<"Cell Pos: "<<hitbar.getPosition().x()<<'\t'<<hitbar.getPosition().y()<<'\t'<<hitbar.getPosition().z()<<endl;
+				cout<<"step pos: "<<steppos.x()<<'\t'<<steppos.y()<<'\t'<<steppos.z()<<endl;
+				cout<<"Relative pos: "<<rpos.x()<<'\t'<<rpos.y()<<'\t'<<rpos.z()<<endl;
+				cout<<"Cell: "<<hitbar.getModule()<<"  "<<hitbar.getDlayer()<<"  "<<hitbar.getSlayer()<<endl;
+			}
+
+			//Get digitalized signal(Q1, Q2, T1, T2) from step
+			//Define: 1 is left, 2 is right, clockwise direction in phi. 
+
+			int sign=-999;
+			if(hitbar.getSlayer()==1) sign = rpos.z()==0 ? 1 : rpos.z()/fabs(rpos.z());
+			else{
+				if(hitbar.getModule()==0 || hitbar.getModule()==1 || hitbar.getModule()==7) sign = rpos.x()==0 ?  1: rpos.x()/fabs(rpos.x());
+				if(hitbar.getModule()==3 || hitbar.getModule()==4 || hitbar.getModule()==5) sign = rpos.x()==0 ? -1:-rpos.x()/fabs(rpos.x());
+				else if(hitbar.getModule()==2) sign = rpos.y()==0 ?  1: rpos.y()/fabs(rpos.y());
+				else if(hitbar.getModule()==6) sign = rpos.y()==0 ? -1:-rpos.y()/fabs(rpos.y());
+			}
+			if(!fabs(sign)) {std::cout<<"ERROR: Wrong bar direction/position!"<<std::endl; continue;}
+
+
+			double Qi_left = en*exp(-(Lbar/2 + sign*sqrt(rpos.Mag2()))/Latt);	
+			double Qi_right = en*exp(-(Lbar/2 - sign*sqrt(rpos.Mag2()))/Latt);
+
+			if(_Debug>=3){
+				cout<<Qi_left<<'\t'<<Qi_right<<endl;
+				cout<<Lbar<<'\t'<<sign*sqrt(rpos.Mag2())<<endl;
+			}
+
+
+			double Ti_left = -1; int looptime=0;
+			while(Ti_left<0){ 
+				// Ti_left = Tinit + rndm.Gaus(nMat*(Lbar/2 + sign*sqrt(rpos.Mag2()))/C, Tres); 
+				Ti_left = Tinit + rndm.Gaus(nMat*(Lbar/2 + sign*sqrt(rpos.Mag2()))/C, Tres) + step_time;  // yyy: add step time 
+				looptime++;
+				if(looptime>500){ std::cout<<"ERROR: Step "<<iCont<<" can not get a positive left-side time!"<<std::endl; break;}
+			}
+			if(looptime>500) continue;		
+			double Ti_right = -1; looptime=0;
+			while(Ti_right<0){ 
+				// Ti_right = Tinit + rndm.Gaus(nMat*(Lbar/2 - sign*sqrt(rpos.Mag2()))/C, Tres); 
+				Ti_right = Tinit + rndm.Gaus(nMat*(Lbar/2 - sign*sqrt(rpos.Mag2()))/C, Tres) + step_time;  // yyy: add step time 
+				looptime++;
+            if(looptime>500){ std::cout<<"ERROR: Step "<<iCont<<" can not get a positive right-side time!"<<std::endl; break;}
+			}
+			if(looptime>500) continue;		
+
+
+			m_step_T1.push_back(Ti_left);
+			m_step_T2.push_back(Ti_right);
+			totQ1 += Qi_left;
+			totQ2 += Qi_right;
+	
+			HitStep stepoutL, stepoutR;
+			stepoutL.setQ(Qi_left); stepoutL.setT(Ti_left);
+			stepoutR.setQ(Qi_right); stepoutR.setT(Ti_right);
+			DigiLvec.push_back(stepoutL);
+			DigiRvec.push_back(stepoutR);
+
+		}
+
+		//Time digitalization
+		//if(_Debug>=2) std::cout<<"Time Digitalize: time at Q >"<<_Qthfrac<<"*totQ"<<std::endl;
+		std::sort(DigiLvec.begin(), DigiLvec.end());
+		std::sort(DigiRvec.begin(), DigiRvec.end());
+		double thQ1=0;
+		double thQ2=0;
+		double thT1, thT2; 
+		for(int iCont=0;iCont<DigiLvec.size();iCont++){
+			thQ1 += DigiLvec[iCont].getQ();
+			if(thQ1>totQ1*_Qthfrac){ 
+				thT1 = DigiLvec[iCont].getT(); 
+				if(_Debug>=3) std::cout<<"Get T1 at index: "<<iCont<<std::endl;
+				break;
+			}
+		}
+		for(int iCont=0;iCont<DigiRvec.size();iCont++){
+			thQ2 += DigiRvec[iCont].getQ();
+			if(thQ2>totQ2*_Qthfrac){ 
+				thT2 = DigiRvec[iCont].getT(); 
+				if(_Debug>=3) std::cout<<"Get T2 at index: "<<iCont<<std::endl;
+				break;
+			}
+		}
+		hitbar.setQ(totQ1, totQ2);
+		hitbar.setT(thT1, thT2);
+    //End bar digitization. 
+
+    //2 hits with double-readout time. 
+    edm4hep::Vector3f m_pos(hitbar.getPosition().X(), hitbar.getPosition().Y(), hitbar.getPosition().Z());
+    auto digiHit1 = caloVec->create();
+    digiHit1.setCellID(hitbar.getcellID());
+    digiHit1.setEnergy(hitbar.getQ1());
+    digiHit1.setTime(hitbar.getT1());
+    digiHit1.setPosition(m_pos);
+    auto digiHit2 = caloVec->create();
+    digiHit2.setCellID(hitbar.getcellID());
+    digiHit2.setEnergy(hitbar.getQ2());
+    digiHit2.setTime(hitbar.getT2());
+    digiHit2.setPosition(m_pos);
+
+    //SimHit - CaloHit association
+    auto rel1 = caloAssoVec->create();
+    rel1.setRec(digiHit1);
+    rel1.setSim(SimHit);
+    rel1.setWeight( hitbar.getQ1()/(hitbar.getQ1()+hitbar.getQ2()) );
+    auto rel2 = caloAssoVec->create();
+    rel2.setRec(digiHit2);
+    rel2.setSim(SimHit);
+    rel2.setWeight( hitbar.getQ2()/(hitbar.getQ1()+hitbar.getQ2()) );
+
+
+    //MCParticle - CaloHit association
+    //float maxMCE = -99.;
+    //edm4hep::MCParticle selMCP; 
+    for(auto iter : MCPEnMap){
+      //if(iter.second>maxMCE){
+      //  maxMCE = iter.second;
+      //  selMCP = iter.first;
+      //}
+      auto rel_MCP1 = caloMCPAssoVec->create();
+      rel_MCP1.setRec(digiHit1);
+      rel_MCP1.setSim(iter.first);
+      rel_MCP1.setWeight(iter.second/digiHit1.getEnergy());
+      auto rel_MCP2 = caloMCPAssoVec->create();
+      rel_MCP2.setRec(digiHit2);
+      rel_MCP2.setSim(iter.first);
+      rel_MCP2.setWeight(iter.second/digiHit2.getEnergy());      
+    }
+
+    //if(selMCP.isAvailable()){
+    //  auto rel_MCP1 = caloMCPAssoVec->create();
+    //  rel_MCP1.setRec(digiHit1);
+    //  rel_MCP1.setSim(selMCP);
+    //  rel_MCP1.setWeight(1.);
+    //  auto rel_MCP2 = caloMCPAssoVec->create();
+    //  rel_MCP2.setRec(digiHit2);
+    //  rel_MCP2.setSim(selMCP);
+    //  rel_MCP2.setWeight(1.);
+    //}
+
+    m_barCol.push_back(hitbar);
+		totE+=(hitbar.getQ1()+hitbar.getQ2())/2;
+		
+    //Temp: write into trees. 
+    if(_writeNtuple){
+      m_simBar_x.push_back(hitbar.getPosition().x());
+      m_simBar_y.push_back(hitbar.getPosition().y());
+      m_simBar_z.push_back(hitbar.getPosition().z());
+      m_simBar_Q1.push_back(hitbar.getQ1());
+      m_simBar_Q2.push_back(hitbar.getQ2());
+      m_simBar_T1.push_back(hitbar.getT1());
+      m_simBar_T2.push_back(hitbar.getT2());
+      m_simBar_module.push_back(hitbar.getModule());
+      m_simBar_stave.push_back(hitbar.getStave());
+      m_simBar_dlayer.push_back(hitbar.getDlayer());
+      m_simBar_part.push_back(hitbar.getPart());
+      m_simBar_slayer.push_back(hitbar.getSlayer());
+      m_simBar_cellID.push_back(hitbar.getcellID());
+    }
+	}
+
+
+  if(_writeNtuple){
+  	t_SimCont->Fill();
+	  t_SimBar->Fill();
+  }
+	if(_Debug>=1) std::cout<<"End Loop: Bar Digitalization!"<<std::endl;
+	std::cout<<"Total Bar Energy: "<<totE<<std::endl;
+
+// yyy_enddigi = clock();
+// double duration_digi = double(yyy_enddigi - yyy_start) / CLOCKS_PER_SEC;
+// // 将时间输出到txt文件中
+// std::ofstream outfile("runtime_ecaldigi.txt", std::ios::app);
+// outfile << _nEvt << "    " << duration_digi << std::endl;
+// outfile.close();
+
+  _nEvt ++ ;
+  //delete SimHitCol, caloVec, caloAssoVec; 
+  m_simhitCol.clear();
+  return StatusCode::SUCCESS;
+}
+
+StatusCode CRDEcalDigiAlg::finalize()
+{
+  if(_writeNtuple){
+  	m_wfile->cd();
+	  t_SimCont->Write();
+  	t_SimBar->Write();
+	  m_wfile->Close();
+    delete m_wfile, t_SimCont, t_SimBar; 
+  }
+
+  info() << "Processed " << _nEvt << " events " << endmsg;
+  delete m_cellIDConverter, m_decoder, m_geosvc;
+  return GaudiAlgorithm::finalize();
+}
+
+StatusCode CRDEcalDigiAlg::MergeHits( const edm4hep::SimCalorimeterHitCollection& m_col, std::vector<edm4hep::SimCalorimeterHit>& m_hits ){
+
+  m_hits.clear(); 
+	std::vector<edm4hep::MutableSimCalorimeterHit> m_mergedhit;
+	m_mergedhit.clear();
+
+	for(int iter=0; iter<m_col.size(); iter++){
+		edm4hep::SimCalorimeterHit m_step = m_col[iter];
+		if(!m_step.isAvailable()){ cout<<"ERROR HIT!"<<endl; continue;}
+		if(m_step.getEnergy()==0 || m_step.contributions_size()<1) continue;
+		unsigned long long cellid = m_step.getCellID();
+		dd4hep::Position hitpos = m_cellIDConverter->position(cellid);
+		edm4hep::Vector3f pos(hitpos.x()*10, hitpos.y()*10, hitpos.z()*10);
+
+		edm4hep::MutableCaloHitContribution conb;
+		conb.setEnergy(m_step.getEnergy());
+		conb.setStepPosition(m_step.getPosition());
+    conb.setParticle( m_step.getContributions(0).getParticle() );
+		conb.setTime(m_step.getContributions(0).getTime());
+
+		edm4hep::MutableSimCalorimeterHit m_hit = find(m_mergedhit, cellid);
+		if(m_hit.getCellID()==0){
+			//m_hit = new edm4hep::SimCalorimeterHit();
+			m_hit.setCellID(cellid);
+			m_hit.setPosition(pos);
+			m_mergedhit.push_back(m_hit);
+		}
+		m_hit.addToContributions(conb);
+		m_hit.setEnergy(m_hit.getEnergy()+m_step.getEnergy());
+	}
+
+  for(auto iter = m_mergedhit.begin(); iter!=m_mergedhit.end(); iter++){
+    edm4hep::SimCalorimeterHit constsimhit = *iter; 
+    m_hits.push_back( constsimhit );  
+  }
+  return StatusCode::SUCCESS; 
+}
+
+
+
+double CRDEcalDigiAlg::GetBarLength(CaloBar& bar){
+	//TODO: reading bar length from geosvc. 
+	if(bar.getSlayer()==1) return 600.;
+	else return 480.-bar.getDlayer()*10.;
+}
+
+/*
+dd4hep::Position CRDEcalDigiAlg::GetCellPos(dd4hep::Position& pos, CaloBar& bar){
+	dd4hep::Position rpos = pos-bar.getPosition();
+	TVector3 vec(0,0,0); 
+	if(bar.getSlayer()==1) vec.SetXYZ(0, 0, floor(rpos.z()/10)*10+5 );
+	else if(bar.getSlayer()==0){
+		if((bar.getModule()==0||bar.getModule()==4) && bar.getDlayer()%2==1) vec.SetXYZ(floor(rpos.x()/10)*10+5,0,0);
+		if((bar.getModule()==0||bar.getModule()==4) && bar.getDlayer()%2==0) vec.SetXYZ(floor((rpos.x()-5)/10)*10+10,0,0);
+		if((bar.getModule()==2||bar.getModule()==6) && bar.getDlayer()%2==1) vec.SetXYZ(0, floor(rpos.y()/10)*10+5,0);
+		if((bar.getModule()==2||bar.getModule()==6) && bar.getDlayer()%2==0) vec.SetXYZ(0, floor((rpos.y()-5)/10)*10+10,0);
+		if(bar.getModule()==1 || bar.getModule()==5){
+			TVector3 unitv(1./sqrt(2), -1./sqrt(2), 0);
+			if(bar.getDlayer()%2==1) vec = (floor(rpos.Dot(unitv)/10)*10+5)*unitv;
+			if(bar.getDlayer()%2==0) vec = (floor((rpos.Dot(unitv)-5)/10)*10+10)*unitv;
+		}
+		if(bar.getModule()==3 || bar.getModule()==7){
+			TVector3 unitv(1./sqrt(2), 1./sqrt(2), 0);
+			if(bar.getDlayer()%2==1) vec = (floor(rpos.Dot(unitv)/10)*10+5)*unitv;
+			if(bar.getDlayer()%2==0) vec = (floor((rpos.Dot(unitv)-5)/10)*10+10)*unitv;
+		}
+	}
+	dd4hep::Position relv(vec.x(), vec.y(), vec.z());
+	return relv+bar.getPosition();
+}
+
+
+edm4hep::MutableSimCalorimeterHit CRDEcalDigiAlg::find(edm4hep::SimCalorimeterHitCollection& m_col, dd4hep::Position& pos){
+   for(int i=0;i<m_col.size();i++){
+    edm4hep::MutableSimCalorimeterHit hit = m_col[i];
+		dd4hep::Position ipos(hit.getPosition().x, hit.getPosition().y, hit.getPosition().z);
+		if(ipos==pos) return hit;
+	}
+   edm4hep::MutableSimCalorimeterHit hit;
+   hit.setCellID(0);
+   return hit;
+}
+*/
+edm4hep::MutableSimCalorimeterHit CRDEcalDigiAlg::find(const std::vector<edm4hep::MutableSimCalorimeterHit>& m_col, unsigned long long& cellid) const{
+   for(int i=0;i<m_col.size();i++){
+		edm4hep::MutableSimCalorimeterHit hit=m_col.at(i);
+		if(hit.getCellID() == cellid) return hit;
+	}
+	edm4hep::MutableSimCalorimeterHit hit ;
+	hit.setCellID(0);
+	return hit;
+}
+
+void CRDEcalDigiAlg::Clear(){
+  totE = -99;
+	m_step_x.clear();
+	m_step_y.clear();
+	m_step_z.clear();
+	m_step_t.clear();   // yyy: clear
+	m_step_E.clear();
+	m_stepBar_x.clear();
+	m_stepBar_y.clear();
+	m_stepBar_z.clear();
+	m_step_T1.clear();
+	m_step_T2.clear();
+	m_simBar_x.clear();
+	m_simBar_y.clear();
+	m_simBar_z.clear();
+	m_simBar_T1.clear();
+	m_simBar_T2.clear();
+	m_simBar_Q1.clear();
+	m_simBar_Q2.clear();
+	m_simBar_module.clear();
+	m_simBar_stave.clear();
+	m_simBar_dlayer.clear();
+	m_simBar_part.clear();
+	m_simBar_slayer.clear();
+  m_simBar_cellID.clear();
+}
+
diff --git a/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.h b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.h
new file mode 100644
index 000000000..cb9be22d6
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.h
@@ -0,0 +1,116 @@
+#ifndef CRD_ECAL_DIGI_ALG_H
+#define CRD_ECAL_DIGI_ALG_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/MutableCaloHitContribution.h"
+#include "edm4hep/MutableSimCalorimeterHit.h"
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoCaloAssociationCollection.h"
+#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"
+#include "CaloBar.h"
+#include "HitStep.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h> 
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h> 
+#include <ctime>
+
+#define C 299.79  // unit: mm/ns
+#define PI 3.141592653
+
+class CRDEcalDigiAlg : public GaudiAlgorithm
+{
+ 
+public:
+ 
+  CRDEcalDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  /** Called at the begin of the job before anything is read.
+   * Use to initialize the processor, e.g. book histograms.
+   */
+  virtual StatusCode initialize() ;
+ 
+  /** Called for every event - the working horse.
+   */
+  virtual StatusCode execute() ; 
+ 
+  /** Called after data processing for clean up.
+   */
+  virtual StatusCode finalize() ;
+
+
+  StatusCode MergeHits(const edm4hep::SimCalorimeterHitCollection& m_col, std::vector<edm4hep::SimCalorimeterHit>& m_hits);  
+
+	double GetBarLength(CaloBar& bar); //TODO: should read from geom file! 
+	edm4hep::MutableSimCalorimeterHit find(const std::vector<edm4hep::MutableSimCalorimeterHit>& m_col, unsigned long long& cellid) const;
+
+	void Clear();
+
+protected:
+
+  SmartIF<IGeomSvc> m_geosvc;
+  //SmartIF<ICRDEcalSvc> m_edmsvc;
+  typedef std::vector<float> FloatVec;
+  typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
+
+	int _nEvt ;
+	float m_length;
+	TRandom3 rndm;
+	TFile* m_wfile;
+	TTree* t_SimCont;
+	TTree* t_SimBar;
+	
+  double totE;
+  FloatVec m_step_t;  // yyy: time of each step
+	FloatVec m_step_x, m_step_y, m_step_z, m_step_E, m_step_T1, m_step_T2, m_stepBar_x, m_stepBar_y, m_stepBar_z;
+	FloatVec m_simBar_x, m_simBar_y, m_simBar_z, m_simBar_T1, m_simBar_T2, m_simBar_Q1, m_simBar_Q2, m_simBar_dlayer, m_simBar_part, m_simBar_stave, m_simBar_slayer, m_simBar_module;
+  std::vector<unsigned long long> m_simBar_cellID;
+
+
+	dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+	dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+	Gaudi::Property<float> m_scale{ this, "Scale", 1 };
+
+  // Input collections
+  DataHandle<edm4hep::SimCalorimeterHitCollection> r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this};
+  mutable Gaudi::Property<std::string> _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"};
+  mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+
+  //Input parameters
+  mutable Gaudi::Property<int>   _writeNtuple{this,  "WriteNtuple", 1, "Write ntuple"};
+  mutable Gaudi::Property<int>   _Nskip{this,  "SkipEvt", 0, "Skip event"};
+  mutable Gaudi::Property<float> _seed{this,   "Seed", 2131, "Random Seed"};
+  mutable Gaudi::Property<int>  _Debug{this,   "Debug", 0, "Debug level"};
+  mutable Gaudi::Property<float> _Eth {this,   "EnergyThreshold", 0.001, "Energy Threshold (/GeV)"};
+  mutable Gaudi::Property<float> r_cali{this,  "CalibrECAL", 1, "Calibration coefficients for ECAL"};
+  mutable Gaudi::Property<float> Latt{this, 	"AttenuationLength", 7000, "Crystal Attenuation Length(mm)"};
+  mutable Gaudi::Property<float> Tres{this, 	"TimeResolution", 0.1, "Crystal time resolution in one side (ns)"};
+  mutable Gaudi::Property<float> nMat{this, 	"MatRefractive", 2.15, "Material refractive index of crystal"};
+  mutable Gaudi::Property<float> Tinit{this, 	"InitalTime", 2, "Start time (ns)"};
+  
+  mutable Gaudi::Property<float> _Qthfrac  {this, 	"ChargeThresholdFrac", 0.05, "Charge threshold fraction"};
+
+
+  // Output collections
+  DataHandle<edm4hep::CalorimeterHitCollection>    w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoCaloAssociationCollection>    w_CaloAssociationCol{"ECALBarrelAssoCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>    w_MCPCaloAssociationCol{"ECALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this};
+};
+
+#endif
diff --git a/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp
new file mode 100644
index 000000000..7491716d3
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp
@@ -0,0 +1,214 @@
+// /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+// // Unit in code: mm, ns. 
+// // NOTE: This digitialization highly matches detector geometry CRDEcalBarrel_v01.
+// // TODO: read geometry info automatically.  
+
+#include "CRDHcalDigiAlg.h" 
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+#include "edm4hep/Cluster.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+
+#include "TVector3.h"
+#include <math.h>
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <map>
+
+// #include <fstream>
+// #include <ctime>
+
+#define C 299.79  // unit: mm/ns
+#define PI 3.141592653
+using namespace std;
+using namespace dd4hep;
+
+DECLARE_COMPONENT( CRDHcalDigiAlg )
+
+CRDHcalDigiAlg::CRDHcalDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+  : GaudiAlgorithm(name, svcLoc),
+    _nEvt(0)
+{
+  
+	// Input collections
+	declareProperty("SimCaloHitCollection", r_SimCaloCol, "Handle of the Input SimCaloHit collection");
+  
+	// Output collections
+	declareProperty("CaloHitCollection", w_DigiCaloCol, "Handle of Digi CaloHit collection");
+	declareProperty("CaloAssociationCollection", w_CaloAssociationCol, "Handle of CaloAssociation collection");
+  declareProperty("CaloMCPAssociationCollection", w_MCPCaloAssociationCol, "Handle of CaloAssociation collection"); 
+}
+
+StatusCode CRDHcalDigiAlg::initialize()
+{
+  if(_writeNtuple){
+    std::string s_outfile = _filename;
+    m_wfile = new TFile(s_outfile.c_str(), "recreate");
+    t_simHit = new TTree("simHit", "simHit");
+    
+    t_simHit->Branch("simHit_x", &m_simHit_x);
+    t_simHit->Branch("simHit_y", &m_simHit_y);
+    t_simHit->Branch("simHit_z", &m_simHit_z);
+    t_simHit->Branch("simHit_E", &m_simHit_E);
+    t_simHit->Branch("simHit_steps", &m_simHit_steps);
+    
+    t_simHit->Branch("simHit_module", &m_simHit_module);
+    t_simHit->Branch("simHit_stave", &m_simHit_stave);
+    t_simHit->Branch("simHit_layer", &m_simHit_layer);
+    t_simHit->Branch("simHit_tower", &m_simHit_tower);
+    t_simHit->Branch("simHit_slice", &m_simHit_slice);
+    t_simHit->Branch("simHit_cellID", &m_simHit_cellID);
+  }
+	std::cout<<"CRDHcalDigiAlg::m_scale="<<m_scale<<std::endl;
+	m_geosvc = service<IGeomSvc>("GeomSvc");
+	if ( !m_geosvc )  throw "CRDHcalDigiAlg :Failed to find GeomSvc ...";
+	dd4hep::Detector* m_dd4hep = m_geosvc->lcdd();
+	if ( !m_dd4hep )  throw "CRDHcalDigiAlg :Failed to get dd4hep::Detector ...";
+	m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep);
+   	m_decoder = m_geosvc->getDecoder(_readoutName);
+	if (!m_decoder) {
+		error() << "Failed to get the decoder. " << endmsg;
+		return StatusCode::FAILURE;
+	}
+
+	//m_edmsvc = service<ICRDEcalSvc>("CRDEcalSvc");
+	//if ( !m_edmsvc )  throw "CRDHcalDigiAlg :Failed to find CRDEcalSvc ...";
+
+	rndm.SetSeed(_seed);
+	std::cout<<"CRDHcalDigiAlg::initialize"<<std::endl;
+	return GaudiAlgorithm::initialize();
+}
+
+StatusCode CRDHcalDigiAlg::execute()
+{
+// clock_t yyy_start, yyy_enddigi;
+// yyy_start = clock(); // 记录开始时间
+
+	if(_nEvt==0) std::cout<<"CRDHcalDigiAlg::execute Start"<<std::endl;
+	std::cout<<"Processing event: "<<_nEvt<<std::endl;
+   	if(_nEvt<_Nskip){ _nEvt++; return StatusCode::SUCCESS; }
+
+	Clear();
+
+ 	const edm4hep::SimCalorimeterHitCollection* SimHitCol =  r_SimCaloCol.get();
+
+	edm4hep::CalorimeterHitCollection* caloVec = w_DigiCaloCol.createAndPut();
+	edm4hep::MCRecoCaloAssociationCollection* caloAssoVec = w_CaloAssociationCol.createAndPut();
+  edm4hep::MCRecoCaloParticleAssociationCollection* caloMCPAssoVec = w_MCPCaloAssociationCol.createAndPut(); 
+
+	if(SimHitCol == 0) 
+	{
+		std::cout<<"not found SimCalorimeterHitCollection"<< std::endl;
+		return StatusCode::SUCCESS;
+	}
+  if(_Debug>=1) std::cout<<"digi, input sim hit size="<< SimHitCol->size() <<std::endl;
+
+
+  for(int isim=0; isim<SimHitCol->size(); isim++){
+
+    auto simhit = SimHitCol->at(isim);
+    if(!simhit.isAvailable()) continue;
+    if(simhit.getEnergy()==0) continue;
+
+    unsigned long long id = simhit.getCellID();
+    double Ehit = simhit.getEnergy();
+    //Energy threshold
+    if(Ehit<_MIPCali*_Eth_Mip) continue;
+
+    //Global calibration. 
+    //TODO: add more digitization terms here. 
+    double Ehit_cali = Ehit*r_cali;
+
+    //Loop contributions to get hit time and MCParticle. 
+    double Thit_ave = 0.;
+    double Ehit_raw = 0.;
+    MCParticleToEnergyWeightMap MCPEnMap; MCPEnMap.clear();
+    for(int iConb=0; iConb<simhit.contributions_size(); ++iConb){
+      auto conb = simhit.getContributions(iConb);
+      if(!conb.isAvailable()) continue;
+      if(conb.getEnergy()==0) continue;
+
+      Thit_ave += conb.getTime();
+      
+      auto mcp = conb.getParticle();
+      MCPEnMap[mcp] += conb.getEnergy();
+      Ehit_raw += conb.getEnergy();
+    }
+    Thit_ave = Thit_ave/simhit.contributions_size();
+    //Create DigiHit
+    auto digiHit = caloVec->create();
+    digiHit.setCellID(id);
+    digiHit.setEnergy(Ehit_cali);
+    digiHit.setTime(Thit_ave);
+    digiHit.setPosition(simhit.getPosition());
+
+    //Create SimHit-DigiHit association. 
+    auto rel = caloAssoVec->create();
+    rel.setRec(digiHit);
+    rel.setSim(simhit);
+    rel.setWeight(1.);
+
+    //Create DigiHit-MCParticle association.
+    for(auto iter : MCPEnMap){
+      auto rel_MC = caloMCPAssoVec->create();
+      rel_MC.setRec(digiHit);
+      rel_MC.setSim(iter.first);
+      rel_MC.setWeight(iter.second/Ehit_raw);
+    }
+
+    if(_writeNtuple){
+      m_simHit_x.push_back(digiHit.getPosition().x);
+      m_simHit_y.push_back(digiHit.getPosition().y);
+      m_simHit_z.push_back(digiHit.getPosition().z);
+      m_simHit_E.push_back(digiHit.getEnergy());
+      m_simHit_steps.push_back(simhit.contributions_size());
+      m_simHit_module.push_back(m_decoder->get(id, "module"));
+      m_simHit_stave.push_back(m_decoder->get(id, "stave"));
+      m_simHit_layer.push_back(m_decoder->get(id, "layer"));
+      m_simHit_slice.push_back(m_decoder->get(id, "slice"));
+      m_simHit_tower.push_back(m_decoder->get(id, "tower"));
+      m_simHit_cellID.push_back(id);
+    }
+  }
+
+	if(_writeNtuple) t_simHit->Fill();
+
+	_nEvt ++ ;
+	return StatusCode::SUCCESS;
+}
+
+StatusCode CRDHcalDigiAlg::finalize()
+{
+  if(_writeNtuple){
+	  m_wfile->cd();
+	  t_simHit->Write();
+    m_wfile->Close();
+	  delete m_wfile, t_simHit; 
+  }
+
+	info() << "Processed " << _nEvt << " events " << endmsg;
+	delete m_cellIDConverter, m_decoder, m_geosvc;
+	return GaudiAlgorithm::finalize();
+}
+
+
+void CRDHcalDigiAlg::Clear(){
+	m_simHit_x.clear();
+	m_simHit_y.clear();
+	m_simHit_z.clear();
+	m_simHit_E.clear();
+	m_simHit_steps.clear();
+	m_simHit_module.clear();
+	m_simHit_stave.clear();
+	m_simHit_layer.clear();
+	m_simHit_slice.clear();
+	m_simHit_tower.clear();
+  m_simHit_cellID.clear();
+}
+
diff --git a/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h
new file mode 100644
index 000000000..6e3571e89
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h
@@ -0,0 +1,100 @@
+#ifndef CRD_HCAL_DIGI_ALG_H 
+#define CRD_HCAL_DIGI_ALG_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/MutableCaloHitContribution.h"
+#include "edm4hep/MutableSimCalorimeterHit.h"
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoCaloAssociationCollection.h"
+#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"
+#include "edm4hep/MCParticle.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h> 
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h> 
+#include <ctime>
+
+#define C 299.79  // unit: mm/ns
+#define PI 3.141592653
+
+class CRDHcalDigiAlg : public GaudiAlgorithm
+{
+ 
+public:
+ 
+  CRDHcalDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  /** Called at the begin of the job before anything is read.
+   * Use to initialize the processor, e.g. book histograms.
+   */
+  virtual StatusCode initialize() ;
+ 
+  /** Called for every event - the working horse.
+   */
+  virtual StatusCode execute() ; 
+ 
+  /** Called after data processing for clean up.
+   */
+  virtual StatusCode finalize() ;
+
+	void Clear();
+
+protected:
+
+  SmartIF<IGeomSvc> m_geosvc;
+  //SmartIF<ICRDEcalSvc> m_edmsvc;
+  typedef std::vector<float> FloatVec;
+  typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
+
+	int _nEvt ;
+	TRandom3 rndm;
+	TFile* m_wfile;
+	TTree* t_simHit;
+
+	FloatVec m_simHit_x, m_simHit_y, m_simHit_z, m_simHit_E, m_simHit_slice, m_simHit_layer, m_simHit_tower, m_simHit_stave, m_simHit_module;
+  std::vector<unsigned long long> m_simHit_cellID;
+  std::vector<int> m_simHit_steps;
+
+
+	dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+	dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+	Gaudi::Property<float> m_scale{ this, "Scale", 1 };
+
+  // Input collections
+  DataHandle<edm4hep::SimCalorimeterHitCollection> r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this};
+  mutable Gaudi::Property<std::string> _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"};
+  mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+
+  //Input parameters
+  mutable Gaudi::Property<int>   _writeNtuple{this,  "WriteNtuple", 1, "Write ntuple"};
+  mutable Gaudi::Property<int>   _Nskip{this,  "SkipEvt", 0, "Skip event"};
+  mutable Gaudi::Property<float> _seed{this,   "Seed", 2131, "Random Seed"};
+  mutable Gaudi::Property<int>  _Debug{this,   "Debug", 0, "Debug level"};
+  mutable Gaudi::Property<float> _MIPCali {this,   "MIPResponse",  0.0005, "MIP response (/GeV)"};
+  mutable Gaudi::Property<float> _Eth_Mip {this,   "MIPThreshold", 0.5, "Energy Threshold (/MIP)"};
+  mutable Gaudi::Property<float> r_cali{this,  "CalibrHCAL", 1, "Global calibration coefficients for HCAL"};
+
+
+  // Output collections
+  DataHandle<edm4hep::CalorimeterHitCollection>    w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoCaloAssociationCollection>    w_CaloAssociationCol{"HCALBarrelAssoCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>    w_MCPCaloAssociationCol{"HCALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this};
+};
+
+#endif
diff --git a/Digitisers/CRDEcalDigi/src/CaloBar.h b/Digitisers/CRDEcalDigi/src/CaloBar.h
new file mode 100644
index 000000000..631925dae
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CaloBar.h
@@ -0,0 +1,56 @@
+#ifndef _CRD_CALOBAR_
+#define _CRD_CALOBAR_
+
+#include <DD4hep/Objects.h>
+#include "TVector3.h"
+
+class CaloBar{
+
+public:
+  CaloBar(unsigned long long _cellID, int _system, int _module, int _slayer, int _dlayer, int _part, int _stave, int _bar, TVector3 _pos, double _Q1, double _Q2, double _T1, double _T2)
+  : cellID(_cellID), system(_system), module(_module), stave(_stave), dlayer(_dlayer), part(_part), slayer(_slayer), bar(_bar), position(_pos), Q1(_Q1), Q2(_Q2), T1(_T1), T2(_T2) {}; 
+  CaloBar() {};
+
+  inline bool operator == (const CaloBar &x) const{
+    return ( (cellID == x.cellID) && getEnergy()==x.getEnergy() );
+  }
+  unsigned long long getcellID() const { return cellID; }
+  int getSystem() const { return system; }
+  int getModule() const { return module; }
+  int getStave()  const { return stave;  }
+  int getDlayer() const { return dlayer; }
+  int getPart()   const { return part;   }
+  int getSlayer() const { return slayer; }
+  int getBar()    const { return bar;    }
+  double getQ1()  const { return Q1;     }
+  double getQ2()  const { return Q2;     }
+  double getT1()  const { return T1;     }
+  double getT2()  const { return T2;     }
+
+  TVector3 getPosition() const { return position; }
+  double getEnergy() const { return (Q1+Q2)/2.; }
+
+  void setcellID(unsigned long long _cellid) { cellID = _cellid; }
+  void setcellID(int _system, int _module, int _stave, int _dlayer, int _part, int _slayer, int _bar) { system=_system; module=_module; stave=_stave; dlayer=_dlayer; part=_part; slayer=_slayer; bar=_bar; }
+  void setPosition( TVector3 posv3) { position.SetXYZ( posv3.x(), posv3.y(), posv3.z() ); }
+  void setQ(double _q1, double _q2) { Q1=_q1; Q2=_q2; }
+  void setT(double _t1, double _t2) { T1=_t1; T2=_t2; }
+
+private:
+	unsigned long long cellID;
+	int system;
+	int module;
+	int stave;
+	int dlayer;
+	int part;
+	int slayer;
+	int bar;
+	TVector3 position;
+	double Q1;      // Q in left readout
+	double Q2;      // Q in right readout;
+	double T1;    // T in left readout;
+	double T2;    // T in right readout;
+
+};
+  
+#endif
diff --git a/Digitisers/CRDEcalDigi/src/CaloHit.h b/Digitisers/CRDEcalDigi/src/CaloHit.h
new file mode 100644
index 000000000..c503f987a
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/CaloHit.h
@@ -0,0 +1,50 @@
+#ifndef _CRD_CALOHIT_
+#define _CRD_CALOHIT_
+
+#include <DD4hep/Objects.h>
+#include "TVector3.h"
+
+class CaloHit{
+
+public:
+  CaloHit(unsigned long long _cellID, int _system, int _module, int _layer, int _stave, int _tower, int _slice, TVector3 _pos, double _E, int _step_numbers)
+  : cellID(_cellID), system(_system), module(_module), layer(_layer), stave(_stave), tower(_tower), slice(_slice), position(_pos), E(_E), step_numbers(_step_numbers) {}; 
+  CaloHit() {};
+
+  inline bool operator == (const CaloHit &x) const{
+    return ( (cellID == x.cellID) && getEnergy()==x.getEnergy() );
+  }
+  unsigned long long getcellID() const { return cellID; }
+  int getSystem() const { return system; }
+  int getModule() const { return module; }
+  int getStave()  const { return stave;  }
+  int getLayer()  const { return layer;  }
+  int getTower()  const { return tower;  }
+  int getSlice()  const { return slice;  }
+
+  TVector3 getPosition() const { return position; }
+  double getEnergy() const { return E; }
+  int getStep_numbers() const { return step_numbers; }
+
+  void setcellID(unsigned long long _cellid) { cellID = _cellid; }
+  void setcellID(int _system, int _module, int _layer, int _stave, int _tower, int _slice) { system=_system; module=_module; stave=_stave; layer=_layer; tower=_tower; slice=_slice;}
+  void setPosition( TVector3 posv3) { position.SetXYZ( posv3.x(), posv3.y(), posv3.z() ); }
+  void setE(double _E) { E=_E; }
+  void setStep_numbers(int _step_numbers) { step_numbers=_step_numbers; }
+
+
+private:
+	unsigned long long cellID;
+	int system;
+	int module;
+	int stave;
+	int layer;
+	int tower;
+	int slice;
+	TVector3 position;
+	double E;
+  int step_numbers;
+
+};
+  
+#endif
diff --git a/Digitisers/CRDEcalDigi/src/HitStep.h b/Digitisers/CRDEcalDigi/src/HitStep.h
new file mode 100644
index 000000000..2ae0715be
--- /dev/null
+++ b/Digitisers/CRDEcalDigi/src/HitStep.h
@@ -0,0 +1,26 @@
+#ifndef HIT_STEP_H
+#define HIT_STEP_H
+
+class HitStep{
+
+public: 
+  HitStep (double _Q, double _T): Q(_Q), T(_T) {};
+  //HitStep (double _Q, double _T) { Q=_Q; T=_T; };
+  HitStep() {};
+
+  void setQ(double _Q) { Q =_Q; }
+  void setT(double _T) { T =_T; }
+
+  double getQ() const { return Q; }
+  double getT() const { return T; }
+  inline bool operator < (const HitStep &x) const { 
+    return T <x.T ;
+  }
+
+private: 
+  double Q;
+  double T;
+
+};
+
+#endif
diff --git a/Examples/options/detsim_FullDet_HCAL.py b/Examples/options/detsim_FullDet_HCAL.py
new file mode 100644
index 000000000..7dfeb66cb
--- /dev/null
+++ b/Examples/options/detsim_FullDet_HCAL.py
@@ -0,0 +1,315 @@
+#!/usr/bin/env python
+
+
+#######  WARNING #########################################
+# 
+# This script can not run in CEPCSW master due to geometry
+# But can be an example for having truth info in track with FullLDCTrackingAlg
+#
+###########################################################
+
+from Gaudi.Configuration import *
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+
+seed = [135]
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = seed
+
+rndmgensvc = RndmGenSvc("RndmGenSvc")
+rndmgensvc.Engine = rndmengine.name()
+
+#geometry_option = "CRD_o1_v01/CRD_o1_v01.xml"
+geometry_option = "CRD_o1_v01/CRD_o1_v01_HCAL.xml"
+#...
+
+if not os.getenv("DETCRDROOT"):
+    print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+    sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+    print("Can't find the compact geometry file: %s"%geometry_path)
+    sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+#ntsvc.Output = ["MyTuples DATAFILE='result.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+
+gun = GtGunTool("GtGunTool")
+gun.Particles = ["gamma","gamma"]
+#gun.Particles = ["nu_e"]
+gun.PositionXs = [0.,0.]
+gun.PositionYs = [0.,0.]
+gun.PositionZs = [0.,0.]
+gun.EnergyMins = [5.,5.] # GeV
+gun.EnergyMaxs = [5.,5.] # GeV
+gun.ThetaMins  = [90.,91.]   # deg
+gun.ThetaMaxs  = [90.,91.]   # deg
+gun.PhiMins    = [0.,0.]   # deg
+gun.PhiMaxs    = [0.,0.]   # deg
+
+
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC240/higgs/exclusive/E240.Pnnh_bb.e0.p0.whizard195/nnh_bb.e0.p0.00001.stdhep"
+
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+# genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+detsimsvc = DetSimSvc("DetSimSvc")
+
+from Configurables import DetSimAlg
+detsimalg = DetSimAlg("DetSimAlg")
+detsimalg.RandomSeeds = seed
+# detsimalg.VisMacs = ["vis.mac"]
+detsimalg.RunCmds = [
+#    "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+    # example_anatool.name()
+#    "ExampleAnaElemTool",
+    "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+detsimalg.PhysicsList = "QGSP_BERT_EMV"
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+#gearsvc.GearXMLFile = "../../Detector/DetCEPCv4/compact/FullDetGear.xml"
+
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import TimeProjectionChamberSensDetTool
+tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
+tpc_sensdettool.TypeOption = 1
+
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+cal_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+cal_sensdettool.CalNamesMergeDisable = ["CaloDetector"]
+cal_sensdettool.CalNamesApplyBirks = ["HcalBarrel"]
+
+# digitization
+vxdhitname  = "VXDTrackerHits"
+sithitname  = "SITTrackerHits"
+tpchitname  = "TPCTrackerHits"
+sethitname  = "SETTrackerHits"
+setspname   = "SETSpacePoints"
+ftdhitname  = "FTDTrackerHits"
+ftdspname   = "FTDSpacePoints"
+from Configurables import PlanarDigiAlg
+digiVXD = PlanarDigiAlg("VXDDigi")
+digiVXD.SimTrackHitCollection = "VXDCollection"
+digiVXD.TrackerHitCollection = vxdhitname
+digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation"
+digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+digiVXD.UsePlanarTag = True
+#digiVXD.OutputLevel = DEBUG
+
+digiSIT = PlanarDigiAlg("SITDigi")
+digiSIT.IsStrip = False
+digiSIT.SimTrackHitCollection = "SITCollection"
+digiSIT.TrackerHitCollection = sithitname
+digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
+digiSIT.ResolutionU = [0.0072]
+digiSIT.ResolutionV = [0.086]
+digiSIT.UsePlanarTag = True
+#digiSIT.OutputLevel = DEBUG
+
+digiSET = PlanarDigiAlg("SETDigi")
+digiSET.IsStrip = False
+digiSET.SimTrackHitCollection = "SETCollection"
+digiSET.TrackerHitCollection = sethitname
+digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
+digiSET.ResolutionU = [0.0072]
+digiSET.ResolutionV = [0.086]
+digiSET.UsePlanarTag = True
+#digiSET.OutputLevel = DEBUG
+
+# two strip tracker hits -> one space point
+from Configurables import SpacePointBuilderAlg
+spSET = SpacePointBuilderAlg("SETBuilder")
+spSET.TrackerHitCollection = sethitname
+spSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
+spSET.SpacePointCollection = setspname
+spSET.SpacePointAssociationCollection = "SETSpacePointAssociation"
+#spSET.OutputLevel = DEBUG
+
+
+digiFTD = PlanarDigiAlg("FTDDigi")
+digiFTD.IsStrip = False
+digiFTD.SimTrackHitCollection = "FTDCollection"
+digiFTD.TrackerHitCollection = ftdhitname
+digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+digiFTD.ResolutionU = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072]
+digiFTD.ResolutionV = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072]
+digiFTD.UsePlanarTag = True
+#digiFTD.OutputLevel = DEBUG
+
+# two strip tracker hits -> one space point
+from Configurables import SpacePointBuilderAlg
+spFTD = SpacePointBuilderAlg("FTDBuilder")
+spFTD.TrackerHitCollection = ftdhitname
+spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+spFTD.SpacePointCollection = ftdspname
+spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation"
+#spFTD.OutputLevel = DEBUG
+
+from Configurables import TPCDigiAlg
+digiTPC = TPCDigiAlg("TPCDigi")
+digiTPC.TPCCollection = "TPCCollection"
+digiTPC.TPCLowPtCollection = "TPCLowPtCollection"
+digiTPC.TPCTrackerHitsCol = tpchitname
+digiTPC.TPCTrackerHitAssCol = "TPCTrackerHitAssociation"
+#digiTPC.OutputLevel = DEBUG
+
+# tracking
+from Configurables import SiliconTrackingAlg
+tracking = SiliconTrackingAlg("SiliconTracking")
+tracking.HeaderCol = "EventHeader"
+tracking.VTXHitCollection = vxdhitname
+tracking.SITHitCollection = sithitname
+tracking.FTDPixelHitCollection = ftdhitname
+tracking.FTDSpacePointCollection = ftdspname
+tracking.SITRawHitCollection = sithitname
+tracking.FTDRawHitCollection = ftdhitname
+tracking.UseSIT = True
+tracking.SmoothOn = False
+tracking.DumpTime = False
+tracking.NDivisionsInTheta = 10
+#tracking.OutputLevel = DEBUG
+
+from Configurables import ForwardTrackingAlg
+forward = ForwardTrackingAlg("ForwardTracking")
+forward.FTDPixelHitCollection = ftdhitname
+forward.FTDSpacePointCollection = ftdspname
+forward.FTDRawHitCollection = ftdhitname
+forward.Chi2ProbCut = 0.0
+forward.HitsPerTrackMin = 3
+forward.BestSubsetFinder = "SubsetSimple"
+forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation",
+                    "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"]
+forward.CriteriaMin = [0,  0.9,  0,  0.995, 0,  0.8, 0,   20,  1.002, 0.1,      0,   0.99, 0,    0.999, 0,   0.99, 0]
+forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08,  99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05]
+forward.DumpTime = False
+#forward.OutputLevel = DEBUG
+
+from Configurables import TrackSubsetAlg
+subset = TrackSubsetAlg("TrackSubset")
+subset.TrackInputCollections = ["ForwardTracks", "SiTracks"]
+subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, ftdspname]
+subset.TrackSubsetCollection = "SubsetTracks"
+subset.DumpTime = False
+#subset.OutputLevel = DEBUG
+
+#TODO: DC reconstruction, as preliminary, use Clupatra like as TPC
+from Configurables import ClupatraAlg
+clupatra = ClupatraAlg("Clupatra")
+clupatra.TPCHitCollection = tpchitname
+#clupatra.DistanceCut = 100.
+#clupatra.MaxDeltaChi2 = 100.
+#clupatra.Chi2Cut = 150.
+#clupatra.OutputLevel = DEBUG
+
+from Configurables import FullLDCTrackingAlg
+full = FullLDCTrackingAlg("FullTracking")
+full.VTXTrackerHits = vxdhitname
+full.SITTrackerHits = sithitname
+full.TPCTrackerHits = tpchitname  # add TPC or DC tracker hit here, if TPC or DC track is set by full.TPCTracks
+full.SETTrackerHits = sethitname
+full.FTDPixelTrackerHits = ftdhitname
+full.FTDSpacePoints = ftdspname
+full.SITRawHits     = sithitname
+full.SETRawHits     = sethitname
+full.FTDRawHits     = ftdhitname
+full.VTXHitRelCol   = "VXDTrackerHitAssociation"
+full.SITHitRelCol   = "SITTrackerHitAssociation"
+full.SETHitRelCol   = "SETSpacePointAssociation"
+full.FTDHitRelCol   = "FTDTrackerHitAssociation"
+full.TPCHitRelCol   = "TPCTrackerHitAssociation"
+full.TPCTracks = "ClupatraTracks" # add standalone TPC or DC track here
+full.SiTracks  = "SubsetTracks"
+full.OutputTracks  = "MarlinTrkTracks"
+full.DumpTime = False
+full.SITHitToTrackDistance = 3.
+full.SETHitToTrackDistance = 5.
+full.MinChi2ProbForSiliconTracks = 0
+#full.OutputLevel = DEBUG
+
+dedxoption = "BetheBlochEquationDedxSimTool"
+from Configurables import DriftChamberSensDetTool
+dc_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+dc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+    dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+    dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+    dedx_simtool.material_Z = 2
+    dedx_simtool.material_A = 4
+    dedx_simtool.scale = 10
+    dedx_simtool.resolution = 0.0001
+
+# output
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "CRD_HCal_GamGam_1deg.root"
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr(
+    TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSET, digiTPC, tracking, forward, subset, full, out],
+    #TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSET, digiTPC, tracking, forward, subset, out],
+    EvtSel = 'NONE',
+    EvtMax = 100,
+    ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc],
+    #ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
+    OutputLevel=INFO
+)
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
index ee5f1feb2..2964b16dc 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -8,14 +8,6 @@
 #include <edm4hep/TrackerHit.h>
 #include <edm4hep/TrackerHit.h>
 #include <edm4hep/Track.h>
-#if __has_include("edm4hep/EDM4hepVersion.h")
-#include "edm4hep/EDM4hepVersion.h"
-#else
-// Copy the necessary parts from  the header above to make whatever we need to work here
-#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
-// v00-09 is the last version without the capitalization change of the track vector members
-#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
-#endif
 
 #include <iostream>
 #include <algorithm>
@@ -124,7 +116,12 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc
   declareProperty("SiTracks", _SiTrackColHdl, "Si Track Collection");
   
   // Input relation collections
-  
+  declareProperty("VTXHitRelCol", _VTXHitRelColHdl, "Vertex TrackerHit association name");  
+  declareProperty("SITHitRelCol", _SITHitRelColHdl, "SIT TrackerHit association name");  
+  declareProperty("SETHitRelCol", _SETHitRelColHdl, "SET TrackerHit association name");  
+  declareProperty("FTDHitRelCol", _FTDHitRelColHdl, "FTD TrackerHit association name");  
+  declareProperty("TPCHitRelCol", _TPCHitRelColHdl, "TPC TrackerHit association name");  
+
   /*
   registerInputCollection(LCIO::LCRELATION,
                           "TPCTracksMCPRelColl",
@@ -140,6 +137,7 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc
   */ 
   // Output track collection
   declareProperty("OutputTracks", _OutputTrackColHdl, "Full LDC track collection name");
+  declareProperty("OutputTrackAssociationCollection", _OutputTrkRelColHdl, "Track - MCParticle association collection name");
 }
 
 
@@ -212,7 +210,8 @@ StatusCode FullLDCTrackingAlg::execute() {
   auto stopwatch = TStopwatch();
   // debug() << endmsg;
   auto outCol = _OutputTrackColHdl.createAndPut();
-  
+  auto outRelCol = _OutputTrkRelColHdl.createAndPut();  
+
   prepareVectors();
   debug() << "************************************PrepareVectors done..." << endmsg;
   
@@ -247,11 +246,20 @@ StatusCode FullLDCTrackingAlg::execute() {
   catch(std::runtime_error& e){
     error() << e.what() << std::endl;
   }
+
+  //Create the association
+  try{
+    CreateTrackAssociation(outCol, outRelCol);
+  }
+  catch(std::runtime_error& e){
+    error() << e.what() << std::endl;
+  }
+
   CleanUp();
   debug() << "Cleanup is done." << endmsg;
   _nEvt++;
   //  getchar();
-  if(m_tuple){
+  if(m_dumpTime && m_tuple){
     m_timeTotal = stopwatch.RealTime()*1000;
     m_tuple->write();
   }
@@ -521,19 +529,11 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
     float z0TrkCand = trkCand->getZ0();
     //    float phi0TrkCand = trkCand->getPhi();
     // FIXME, fucd
-#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
     int nhits_in_vxd = track.getSubdetectorHitNumbers(0);
     int nhits_in_ftd = track.getSubdetectorHitNumbers(1);
     int nhits_in_sit = track.getSubdetectorHitNumbers(2);
     int nhits_in_tpc = track.getSubdetectorHitNumbers(3);
     int nhits_in_set = track.getSubdetectorHitNumbers(4);
-#else
-    int nhits_in_vxd = track.getSubDetectorHitNumbers(0);
-    int nhits_in_ftd = track.getSubDetectorHitNumbers(1);
-    int nhits_in_sit = track.getSubDetectorHitNumbers(2);
-    int nhits_in_tpc = track.getSubDetectorHitNumbers(3);
-    int nhits_in_set = track.getSubDetectorHitNumbers(4);
-#endif
     //int nhits_in_vxd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - 2 ];
     //int nhits_in_ftd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 2 ];
     //int nhits_in_sit = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SIT - 2 ];
@@ -595,7 +595,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
       debug() << " Add Track to final Collection: ID = " << track.id() << " for trkCand "<< trkCand << endmsg;
     }
   }
-  if(m_tuple) m_timeKalman = stopwatch.RealTime()*1000;
+  if(m_dumpTime && m_tuple) m_timeKalman = stopwatch.RealTime()*1000;
   // streamlog_out(DEBUG5) << endmsg;
   debug() << "Number of accepted " << _OutputTrackColHdl.fullKey() << " = " << nTotTracks << endmsg;
   debug() << "Total 4-momentum of " << _OutputTrackColHdl.fullKey() << " : E = " << eTot
@@ -1166,18 +1166,13 @@ void FullLDCTrackingAlg::prepareVectors() {
       trackExt->setNDF(tpcTrack.getNdf());
       trackExt->setChi2(tpcTrack.getChi2());
       for (int iHit=0;iHit<nHits;++iHit) {
-	edm4hep::TrackerHit hit = tpcTrack.getTrackerHits(iHit);
-	if (!hit.isAvailable()) {
-	  error() << "Tracker hit not available" << endmsg;
-	  continue;
-	}
+	edm4hep::TrackerHit hit = tpcTrack.getTrackerHits(iHit);//hitVec[iHit];
+	if(!hit.isAvailable()) error() << "Tracker hit not available" << endmsg;
 	//info() << "hit " << hit.id() << " " << hit.getCellID() << " " << hit.getPosition()[0] << " " << hit.getPosition()[1] << " " << hit.getPosition()[2] << endmsg;
 	auto it = mapTrackerHits.find(hit);
-	if (it==mapTrackerHits.end()) {
-	  error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
-	  continue;
-	}
-        TrackerHitExtended* hitExt = it->second;
+	if(it==mapTrackerHits.end()) error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
+	else continue;
+        TrackerHitExtended * hitExt = it->second;
 	//info() << hit.id() << " " << hitExt << endmsg;
         hitExt->setTrackExtended( trackExt );
         trackExt->addTrackerHitExtended( hitExt );
@@ -1236,17 +1231,8 @@ void FullLDCTrackingAlg::prepareVectors() {
       char strg[200];
       HelixClass helixSi;
       for (int iHit=0;iHit<nHits;++iHit) {
-	edm4hep::TrackerHit hit = siTrack.getTrackerHits(iHit);
-	if (!hit.isAvailable()) {
-	  error() << "Tracker hit not available" << endmsg;
-	  continue;
-	}
-	auto it = mapTrackerHits.find(hit);
-        if (it==mapTrackerHits.end()) {
-	  error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
-	  continue;
-	}
-        TrackerHitExtended* hitExt = it->second;
+	edm4hep::TrackerHit hit = siTrack.getTrackerHits(iHit);//hitVec[iHit];
+        TrackerHitExtended * hitExt = mapTrackerHits[hit];
         hitExt->setTrackExtended( trackExt );
         
         trackExt->addTrackerHitExtended( hitExt );
@@ -1548,8 +1534,8 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
   int nTPCHits = int(tpcHitVec.size());
   int nHits = nTPCHits + nSiHits;
   
-  //debug() << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg;
-  //debug() << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg;
+  //std::cout << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg;
+  //std::cout << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg;
   
   TrackerHitVec trkHits;
   trkHits.reserve(nHits);
@@ -1762,7 +1748,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
       tpcHitInFit.push_back(tpcHitVec[i]);
     }
   }
-
+  
   debug() << "FullLDCTrackingAlg::CombineTracks: Check for Silicon Hit rejections ... " << endmsg;
   
   if ( (int)siOutliers.size() > _maxAllowedSiHitRejectionsForTrackCombination ) {
@@ -3819,7 +3805,7 @@ HelixClass * FullLDCTrackingAlg::GetExtrapolationHelix( TrackExtended * track) {
           ts_at_calo = getTrackStateAt(trk_lcio, 4/*lcio::TrackState::AtCalorimeter*/);
 
           debug() << "FullLDCTrackingAlg::GetExtrapolationHelix set ts_at_calo with ref_z = " << z_ref << endmsg;
-	}
+	      }
       }
     }
   }
@@ -5013,3 +4999,88 @@ void FullLDCTrackingAlg::checkTrackState(int type){
     }
   }
 }
+
+void FullLDCTrackingAlg::CreateTrackAssociation(edm4hep::TrackCollection* trkCol, edm4hep::MCRecoTrackParticleAssociationCollection* relCol){
+
+  const edm4hep::MCRecoTrackerAssociationCollection* VTXHitRelCol = _VTXHitRelColHdl.get(); 
+  const edm4hep::MCRecoTrackerAssociationCollection* SITHitRelCol = _SITHitRelColHdl.get(); 
+  const edm4hep::MCRecoTrackerAssociationCollection* SETHitRelCol = _SETHitRelColHdl.get(); 
+  const edm4hep::MCRecoTrackerAssociationCollection* FTDHitRelCol = _FTDHitRelColHdl.get(); 
+  const edm4hep::MCRecoTrackerAssociationCollection* TPCHitRelCol = _TPCHitRelColHdl.get(); 
+
+
+  //Loop for the tracks
+  for(int itrk=0; itrk<trkCol->size(); itrk++){
+    edm4hep::Track iTrk = trkCol->at(itrk);
+    if(!iTrk.isAvailable()) continue;
+
+    std::map<const edm4hep::MCParticle, int> MCParticleToNhitMap;
+    //Loop for trackerhits in the track
+    int nHit = iTrk.trackerHits_size();
+    for(int iHit=0; iHit<nHit; ++iHit){
+      const edm4hep::TrackerHit hit = iTrk.getTrackerHits(iHit);
+      if(!hit.isAvailable()) continue;
+
+      //Get SimTrackerHit from Association: loop all kinds of associations
+      bool findSimHit = false; 
+      edm4hep::SimTrackerHit p_simHit;
+      for(int irel = 0; irel<VTXHitRelCol->size() && !findSimHit; ++irel){
+        if(hit == VTXHitRelCol->at(irel).getRec()){
+          p_simHit = VTXHitRelCol->at(irel).getSim();
+          findSimHit = true;
+          break;
+        }
+      }
+
+      for(int irel = 0; irel<SITHitRelCol->size() && !findSimHit; ++irel){
+        if(hit == SITHitRelCol->at(irel).getRec()){
+          p_simHit = SITHitRelCol->at(irel).getSim();
+          findSimHit = true;
+          break;
+        } 
+      }      
+
+      for(int irel = 0; irel<SETHitRelCol->size() && !findSimHit; ++irel){
+        if(hit == SETHitRelCol->at(irel).getRec()){
+          p_simHit = SETHitRelCol->at(irel).getSim();
+          findSimHit = true;
+          break;
+        }
+      }
+
+      for(int irel = 0; irel<FTDHitRelCol->size() && !findSimHit; ++irel){
+        if(hit == FTDHitRelCol->at(irel).getRec()){
+          p_simHit = FTDHitRelCol->at(irel).getSim();
+          findSimHit = true;
+          break;
+        }
+      }
+
+      for(int irel = 0; irel<TPCHitRelCol->size() && !findSimHit; ++irel){
+        if(hit == TPCHitRelCol->at(irel).getRec()){
+          p_simHit = TPCHitRelCol->at(irel).getSim();
+          findSimHit = true;
+          break;
+        }
+      }
+
+      //Get MCP from SimTrackerHit
+      if(findSimHit && p_simHit.isAvailable()){
+        const edm4hep::MCParticle mcp = p_simHit.getMCParticle();
+        if(MCParticleToNhitMap.find(mcp)==MCParticleToNhitMap.end()) MCParticleToNhitMap[mcp]=1;
+        else MCParticleToNhitMap[mcp]++;
+      }
+    } //End loop tracker hits. 
+
+    for(auto& iter : MCParticleToNhitMap){
+      auto rel = relCol->create();
+      rel.setRec(iTrk);
+      rel.setSim(iter.first);
+      rel.setWeight( (double)iter.second/nHit  );
+    }
+
+  } //End loop tracks. 
+}
+
+
+
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
index cfc804879..e7d6c8791 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
@@ -23,6 +23,8 @@
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/TrackerHitCollection.h"
 #include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
 
 #include <UTIL/BitField64.h>
 #include <UTIL/ILDConf.h>
@@ -297,6 +299,8 @@ class FullLDCTrackingAlg : public GaudiAlgorithm {
   void AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, float dcut, int refit);
   
   void CreateExtrapolations();
+
+  void CreateTrackAssociation(edm4hep::TrackCollection* trkCol, edm4hep::MCRecoTrackParticleAssociationCollection* relCol);
   
   void CleanUpExtrapolations();
   
@@ -511,10 +515,17 @@ class FullLDCTrackingAlg : public GaudiAlgorithm {
   DataHandle<edm4hep::TrackerHitCollection> _inFTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this};
   //DataHandle<edm4hep::SimTrackerHitCollection> _inVTXRawColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this};
 
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _VTXHitRelColHdl{"VTXHitRelCol", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _SITHitRelColHdl{"SITHitRelCol", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _SETHitRelColHdl{"SETHitRelCol", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _FTDHitRelColHdl{"FTDHitRelCol", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _TPCHitRelColHdl{"TPCHitRelCol", Gaudi::DataHandle::Reader, this};
+
   DataHandle<edm4hep::TrackCollection> _TPCTrackColHdl{"ClupatraTracks", Gaudi::DataHandle::Reader, this};
   DataHandle<edm4hep::TrackCollection> _SiTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this};
 
   DataHandle<edm4hep::TrackCollection> _OutputTrackColHdl{"MarlinTrkTracks", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> _OutputTrkRelColHdl{"MarlinTrkAssociation", Gaudi::DataHandle::Writer, this};
 
   NTuple::Tuple*       m_tuple;
   NTuple::Item<float>  m_timeTotal;