diff --git a/Digitization/DigiCalo/src/EcalDigiAlg.cpp b/Digitization/DigiCalo/src/EcalDigiAlg.cpp
index 7e0af0cb28602c356206e7a6b5ce5a24504493d3..53e5ecf2d51aa32f24859256bca1a17daaaac259 100644
--- a/Digitization/DigiCalo/src/EcalDigiAlg.cpp
+++ b/Digitization/DigiCalo/src/EcalDigiAlg.cpp
@@ -321,6 +321,9 @@ StatusCode EcalDigiAlg::execute()
       		totQ1_Truth += en;
       
       		auto mcp = conb.getParticle();
+          if( mcp.parents_size()>0 && (fabs(mcp.getVertex().z)>EcalHalfZ || sqrt(mcp.getVertex().x*mcp.getVertex().x + mcp.getVertex().y*mcp.getVertex().y)>EcalInnerR) ){
+            mcp = mcp.getParents(0);
+          }
       		MCPEnMap[mcp] += en;
       		TVector3 steppos(conb.getStepPosition().x, conb.getStepPosition().y, conb.getStepPosition().z);
       		TVector3 rpos = steppos-hitbar.getPosition();
@@ -576,8 +579,17 @@ StatusCode EcalDigiAlg::execute()
       	rel_MCP2.setSim(iter.first);
       	rel_MCP2.setWeight(iter.second/SimHit.getEnergy());      
       	}
-      
-      
+/*
+cout<<"  In SimHit #"<<i<<": linked truth particle "<<MCPEnMap.size()<<": "<<endl;
+for(auto iter : MCPEnMap){
+  cout<<"    Particle "<<iter.first.getPDG()<<", gStatus "<<iter.first.getGeneratorStatus()<<", sStatus "<<iter.first.getSimulatorStatus();
+  cout<<", vertex ("<<iter.first.getVertex().x<<", "<<iter.first.getVertex().y<<", "<<iter.first.getVertex().z<<") ";
+  cout<<", parent size "<<iter.first.parents_size();
+  if(iter.first.parents_size()>0) cout<<" parent ID "<<iter.first.getParents(0).getPDG();
+  cout<<", weight "<<iter.second/SimHit.getEnergy()<<endl;
+  cout<<endl;
+}
+*/
       	//if(selMCP.isAvailable()){
       	//  auto rel_MCP1 = caloMCPAssoVec->create();
       	//  rel_MCP1.setRec(digiHit1);
diff --git a/Digitization/DigiCalo/src/EcalDigiAlg.h b/Digitization/DigiCalo/src/EcalDigiAlg.h
index 4a7a8d13aeb34ef442c529b565b7ffdbdfd47b1f..894f0d02a2d146e80b04932f58679197c7badd64 100644
--- a/Digitization/DigiCalo/src/EcalDigiAlg.h
+++ b/Digitization/DigiCalo/src/EcalDigiAlg.h
@@ -74,6 +74,9 @@ protected:
   typedef std::vector<float> FloatVec;
   typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
 
+  const float EcalInnerR = 1830.;
+  const float EcalHalfZ = 2900.;
+
 	int _nEvt ;
 	float m_length;
 	TRandom3 rndm;
diff --git a/Reconstruction/RecPFACyber/include/CyberDataCol.h b/Reconstruction/RecPFACyber/include/CyberDataCol.h
index 1d6854f6239eaa240bb5172984d17679f648b748..161d824ee88e2d303d7c78dc0254620e33a692fa 100644
--- a/Reconstruction/RecPFACyber/include/CyberDataCol.h
+++ b/Reconstruction/RecPFACyber/include/CyberDataCol.h
@@ -9,8 +9,8 @@
 //  Contact: guofangyi@ihep.ac.cn,
 //           sunss@ihep.ac.cn
 //=============================================================
-#ifndef _PANDORAPLUS_DATA_H
-#define _PANDORAPLUS_DATA_H
+#ifndef _CYBER_DATA_H
+#define _CYBER_DATA_H
 #include <iostream>
 #include <algorithm>
 #include <map>
@@ -32,6 +32,7 @@
 #include "edm4hep/MCParticle.h"
 #include "edm4hep/Track.h"
 #include "edm4hep/TrackCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
 #include "edm4hep/CalorimeterHit.h"
 #include "edm4hep/CalorimeterHitCollection.h"
 #include "edm4hep/Vertex.h"
@@ -45,9 +46,10 @@
 #include "edm4hep/ReconstructedParticleCollection.h"
 #include "edm4hep/MCRecoCaloAssociation.h"
 #include "edm4hep/MCRecoTrackerAssociation.h"
-#include "edm4hep/MCRecoParticleAssociationCollection.h"
-#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"
-#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/MCRecoParticleAssociationCollection.h"         // PFO - MCP
+#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"     // Hit - MCP 
+#include "edm4hep/MCRecoClusterParticleAssociationCollection.h"  // Cluster - MCP
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"    // Track - MCP
 
 #define PI 3.141592653
 //#define C 299.79  // unit: mm/ns
diff --git a/Reconstruction/RecPFACyber/include/CyberPFAlg.h b/Reconstruction/RecPFACyber/include/CyberPFAlg.h
index 6f1d195ee7fd38f79ca2e882d4d640249e881638..a3bb35e8fe09f2b205a0d0118f64e571857aa4be 100644
--- a/Reconstruction/RecPFACyber/include/CyberPFAlg.h
+++ b/Reconstruction/RecPFACyber/include/CyberPFAlg.h
@@ -7,8 +7,8 @@
 //  Contact: guofangyi@ihep.ac.cn,
 //           sunss@ihep.ac.cn
 //=============================================================
-#ifndef PANDORAPLUS_ALG_H
-#define PANDORAPLUS_ALG_H
+#ifndef CYBERPFA_ALG_H
+#define CYBERPFA_ALG_H
 
 #include <string>
 #include "k4FWCore/DataHandle.h"
@@ -136,6 +136,7 @@ protected:
 
   //---Readin collections
   typedef DataHandle<edm4hep::TrackCollection>                          TrackType; 
+  //typedef DataHandle<edm4hep::SimCalorimeterHitCollection>              SimCaloType; 
   typedef DataHandle<edm4hep::CalorimeterHitCollection>                 CaloType; 
   typedef DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>  CaloParticleAssoType; 
   DataHandle<edm4hep::MCParticleCollection>* r_MCParticleCol; 
@@ -143,6 +144,7 @@ protected:
   std::vector<TrackType*> r_TrackCols; 
   //std::vector<CaloType*>  r_ECalHitCols; 
   //std::vector<CaloType*>  r_HCalHitCols; 
+  //std::vector<SimCaloType*>  r_SimCaloHitCols; 
   std::vector<CaloType*>  r_CaloHitCols; 
   std::map<std::string, CaloParticleAssoType*> map_CaloMCPAssoCols;
   DataHandle<edm4hep::RecTofCollection>* r_TofCol;
@@ -181,9 +183,13 @@ protected:
   //Gaudi::Property< std::string > name_EcalCore   {this, "OutputEcalCore", "EcalCore"};
   //Gaudi::Property< std::string > name_HcalCluster{this, "OutputHcalCluster", "HcalCluster"};
   Gaudi::Property< std::string > name_PFObject   {this, "OutputPFO", "outputPFO"};
-  DataHandle<edm4hep::ReconstructedParticleCollection>      w_ReconstructedParticleCollection {"PandoraPFOs"    ,Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::ReconstructedParticleCollection>      w_ReconstructedParticleCollection {"CyberPFO"    ,Gaudi::DataHandle::Writer, this};
   std::map<std::string, DataHandle<edm4hep::ClusterCollection>* > w_ClusterCollection;
 
+  DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>     w_MCPHitAssoCol{"RecCaloHitParticleAssoCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection>  w_MCPClusterAssoCol{"RecCaloClusterAssoCol", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoParticleAssociationCollection>         w_MCPRecoPFOAssoCol{"RecPFOAssoCol", Gaudi::DataHandle::Writer, this};
+
 
   //For Ana
   Gaudi::Property<bool>  m_WriteAna {this, "WriteAna", false, "Write Ntuples for analysis"};
diff --git a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
index 5d4eeb8dfca89e1b9c56e9e7f66ad16721db2e5a..59d17c0133bd669cd516c39ba93da404b146148f 100644
--- a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
+++ b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
@@ -88,6 +88,7 @@ namespace Cyber{
     static int NbarZ;
     static float barsize;
     static float ecal_innerR;
+    static float ecal_innerZ;
     static float ecal_endcap_deadarea;
     static float ecal_endcap_barsize;
 
diff --git a/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h b/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h
index 401e7855956ce04c63ec3cc3bd2deac53a222861..c771296b7bd4d9bc4d3a2f10ce7d82a618bbe830 100644
--- a/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h
+++ b/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h
@@ -16,13 +16,17 @@ namespace Cyber{
     OutputCreator( const Settings& m_settings);
     ~OutputCreator() {};
 
-    StatusCode CreateOutputCollections( CyberDataCol& m_DataCol, 
-                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler, 
-                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler, 
-                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler, 
-                                        DataHandle<edm4hep::TrackCollection>& m_outTrkHandler, 
-                                        std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler, 
-                                        DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler );
+    StatusCode CreateOutputCollections( CyberDataCol& m_DataCol,
+                                        std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol,
+                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler,
+                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler,
+                                        DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler,
+                                        DataHandle<edm4hep::TrackCollection>& m_outTrkHandler,
+                                        std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler,
+                                        DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler,
+                                        DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>& m_outMCPHitLinkHandler,
+                                        DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection>& m_outMCPClusterLinkHandler,
+                                        DataHandle<edm4hep::MCRecoParticleAssociationCollection>& m_outMCPPFOLinkHandler );
 
 
     StatusCode Reset() { return StatusCode::SUCCESS; }
@@ -31,7 +35,7 @@ namespace Cyber{
 
   private: 
     const Cyber::Settings   settings;
-  
+ 
   };
 };
 #endif
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
index e7147853d4521b0c4cbf38486a26ce6cf511812b..bde383c656de0282ffbc38057f2f92c55a08beca 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
@@ -95,6 +95,7 @@ StatusCode HcalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
   //cout<<"  Isolated hit size: "<<Nisohits<<endl;
  
 //   m_datacol.map_CaloCluster[settings.map_stringPars["OutputCluster"]] = m_clusterCol;
+  for(int ic=0; ic<m_clusterCol.size(); ic++) m_clusterCol[ic]->getLinkedMCPfromHit();
   m_datacol.map_CaloCluster[settings.map_stringPars["OutputHCALClusters"]]= m_clusterCol;
   return StatusCode::SUCCESS;
 };
diff --git a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
index f8386eab63e3b44fe11d4e9d354a064795cdaac3..3cc18c0b6ba88e96c1a7a1f2f0908dff8c27608d 100644
--- a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
+++ b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
@@ -9,8 +9,8 @@
 //  Contact: guofangyi@ihep.ac.cn,
 //           sunss@ihep.ac.cn
 //=============================================================
-#ifndef _PANDORAPLUS_DATA_C
-#define _PANDORAPLUS_DATA_C
+#ifndef _CYBER_DATA_C
+#define _CYBER_DATA_C
 
 #include "CyberDataCol.h"
 
diff --git a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
index 42f5421ca7a296badc5b9922f4b1f5ff0886babc..b3cbcfc1a8533731f8920f3ac5bd7c1b3258a23d 100644
--- a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
@@ -8,8 +8,8 @@
 //           sunss@ihep.ac.cn
 //=============================================================
 
-#ifndef PANDORAPLUS_ALG_C
-#define PANDORAPLUS_ALG_C
+#ifndef CYBERPFA_ALG_C
+#define CYBERPFA_ALG_C
 
 #include "CyberPFAlg.h"
 
@@ -31,6 +31,7 @@ int Cyber::CaloUnit::NbarZ = 24;
 //int Cyber::CaloUnit::over_module_set = 2;
 float Cyber::CaloUnit::barsize = 15.2; //mm
 float Cyber::CaloUnit::ecal_innerR = 1830;  //mm
+float Cyber::CaloUnit::ecal_innerZ = 2900;  //mm
 float Cyber::CaloUnit::ecal_endcap_deadarea = 10.5; //mm, a bit larger than real value 8.5 mm in geometry
 float Cyber::CaloUnit::ecal_endcap_barsize = 16.2; //mm, a bit larger than real value 15.2 mm in geometry
 
@@ -48,7 +49,9 @@ CyberPFAlg::CyberPFAlg(const std::string& name, ISvcLocator* svcLoc)
   declareProperty("RecTrackCollection", w_RecTrkCol, "Handle of Reconstructed Tracks linked to PFO"); 
   declareProperty("RecoPFOCollection", w_ReconstructedParticleCollection, "Handle of Reconstructed PFO collection");   
   // declareProperty("RecoVtxCollection", w_VertexCollection, "Handle of Reconstructed vertex collection");   
-  // declareProperty("MCRecoPFOAssociationCollection", w_MCRecoParticleAssociationCollection, "Handle of MC-RecoPFO association collection");   
+ declareProperty("MCRecoHitAssociationCollection", w_MCPHitAssoCol, "Handle of MC-CaloHit association collection");   
+ declareProperty("MCRecoClusterAssociationCollection", w_MCPClusterAssoCol, "Handle of MC-CaloCluster association collection");   
+ declareProperty("MCRecoPFOAssociationCollection", w_MCPRecoPFOAssoCol, "Handle of MC-RecoPFO association collection");   
 
 }
 
@@ -90,6 +93,7 @@ StatusCode CyberPFAlg::initialize()
   m_CaloHitsCreatorSettings.map_stringVecPars["CaloHitCollections"] = name_CaloHits;
   m_CaloHitsCreatorSettings.map_stringPars["EcalType"] = m_EcalType.value();
 
+  m_OutputCreatorSettings.map_stringVecPars["CaloHitCollections"] = name_CaloHits;
   m_OutputCreatorSettings.map_stringPars["OutputPFO"] = name_PFObject.value();
   m_OutputCreatorSettings.map_boolPars["UseTruthTrk"] = m_useMCPTrk.value();
   m_OutputCreatorSettings.map_boolPars["UseTruthMatchTrk"] = m_useTruthMatchTrk.value();
@@ -98,6 +102,9 @@ StatusCode CyberPFAlg::initialize()
   m_OutputCreatorSettings.map_floatPars["HCALChargedCalib"] = m_HcalChargedCalib.value();
   m_OutputCreatorSettings.map_floatPars["ECALNeutralCalib"] = m_EcalNeutralCalib.value();
   m_OutputCreatorSettings.map_floatPars["HCALNeutralCalib"] = m_HcalNeutralCalib.value();
+  m_OutputCreatorSettings.map_floatPars["EcalHalfZ"] = Cyber::CaloUnit::ecal_innerZ;
+  m_OutputCreatorSettings.map_floatPars["EcalInnerR"] = Cyber::CaloUnit::ecal_innerR;
+  m_OutputCreatorSettings.map_floatPars["HitSize"] = Cyber::CaloUnit::barsize;
 
   //Initialize Creators
   m_pMCParticleCreator = new MCParticleCreator( m_pMCParticleCreatorSettings );
@@ -117,6 +124,17 @@ StatusCode CyberPFAlg::initialize()
 
 
   //---Calo Hits---
+  //------Sim hits ------
+  //for(auto& _ecal : name_SimEcalHits){
+  //  if(!_ecal.empty()){
+  //    r_SimCaloHitCols.push_back( new SimCaloType(_ecal, Gaudi::DataHandle::Reader, this) );
+  //}}
+  //for(auto& _hcal : name_SimHcalHits){
+  //  if(!_hcal.empty()){
+  //    r_SimCaloHitCols.push_back( new SimCaloType(_hcal, Gaudi::DataHandle::Reader, this) );
+  //}}  
+
+  //------Digi hits ------
   for(auto& _ecal : name_EcalHits){
     if(!_ecal.empty()){ 
       //r_ECalHitCols.push_back( new CaloType(_ecal, Gaudi::DataHandle::Reader, this) );
@@ -748,12 +766,17 @@ StatusCode CyberPFAlg::execute()
   m_algorithmManager.RunAlgorithm( m_DataCol );
 
   m_pOutputCreator->CreateOutputCollections( m_DataCol, 
+                                             map_CaloMCPAssoCols,
                                              w_RecEcalCol, 
                                              w_RecCoreCol, 
                                              w_RecHcalCol,
                                              w_RecTrkCol, 
                                              w_ClusterCollection, 
-                                             w_ReconstructedParticleCollection);
+                                             w_ReconstructedParticleCollection, 
+                                             w_MCPHitAssoCol, 
+                                             w_MCPClusterAssoCol,
+                                             w_MCPRecoPFOAssoCol );
+
 
 // yyy_endrec = clock();  // 閲嶅缓缁撴潫鐨勬椂闂�
 
@@ -1634,7 +1657,7 @@ StatusCode CyberPFAlg::execute()
     for(int ip=0; ip<m_DataCol.map_PFObjects["outputPFO"].size(); ip++)
       m_pfobjects.push_back(m_DataCol.map_PFObjects["outputPFO"][ip].get());
    
-   
+//cout<<"PFO size: "<<m_pfobjects.size()<<endl;
     for(int ip=0; ip<m_pfobjects.size(); ip++){
       std::vector<const Track*> t_tracks = m_pfobjects[ip]->getTracks();
       std::vector<const Calo3DCluster*> t_ecal_clusters = m_pfobjects[ip]->getECALClusters();
@@ -1658,7 +1681,7 @@ StatusCode CyberPFAlg::execute()
           pfo_trk_location.push_back( AllTrackStates[istate].location );
         }
       }
-   
+//cout<<"  In PFO #"<<ip<<": Ecal cluster size "<<t_ecal_clusters.size()<<endl;   
       for(int ie=0; ie<t_ecal_clusters.size(); ie++){
         pfo_ecal_tag.push_back(ip);
         pfo_ecal_clus_x.push_back(t_ecal_clusters[ie]->getShowerCenter().x());
@@ -1670,7 +1693,17 @@ StatusCode CyberPFAlg::execute()
         //if (tmp_phi < 0) tmp_phi += 360.0;
         //double tmp_theta = std::atan2(t_ecal_clusters[ie]->getShowerCenter().z(), t_ecal_clusters[ie]->getShowerCenter().Perp())* 180.0 / M_PI + 90; 
         pfo_ecal_clus_Escale.push_back(t_ecal_clusters[ie]->getLongiE());
-   
+/*
+double tmp_totE = 0;
+cout<<"    Cluster #"<<ie<<": energy "<<t_ecal_clusters[ie]->getLongiE()<<", all hitE "<<t_ecal_clusters[ie]->getHitsE()<<", hit size "<<t_ecal_clusters[ie]->getCaloHits().size()<<endl;
+for(int ahit=0; ahit<t_ecal_clusters[ie]->getCaloHits().size(); ahit++){
+  const Cyber::CaloHit* p_hit = t_ecal_clusters[ie]->getCaloHits()[ahit];
+  printf("      Hit #%d: position (%.3f, %.3f, %.3f), En %.3f \n", ahit, p_hit->getPosition().x(), p_hit->getPosition().y(), p_hit->getPosition().z(), p_hit->getEnergy());  
+  tmp_totE += p_hit->getEnergy();
+}
+cout<<"    tot Hit E "<<tmp_totE*t_ecal_clusters[ie]->getEnergyScale()<<endl;
+cout<<endl;
+*/
       }
       for(int ih=0; ih<t_hcal_clusters.size(); ih++){
         pfo_hcal_tag.push_back(ip);
@@ -1723,15 +1756,13 @@ StatusCode CyberPFAlg::finalize()
 
   delete r_MCParticleCol;
   for(auto iter : r_TrackCols) delete iter; 
-  //for(auto iter : r_ECalHitCols) delete iter; 
-  //for(auto iter : r_HCalHitCols) delete iter; 
+  //for(auto iter : r_SimCaloHitCols) delete iter; 
   for(auto iter : r_CaloHitCols) delete iter; 
   for(auto iter : map_CaloMCPAssoCols) delete iter.second;
   r_TrackCols.clear();
   for(auto iter : w_ClusterCollection) delete iter.second;
   w_ClusterCollection.clear(); 
-  //r_ECalHitCols.clear();
-  //r_HCalHitCols.clear();
+  //r_SimCaloHitCols.clear();
   r_CaloHitCols.clear();
   //m_energycorsvc->finalize();
   delete m_cellIDConverter, m_geosvc;
diff --git a/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc b/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc
index eb642f53d5efb7f7ffd1e246f93ff52b4b8fc747..4bbe660eb17a98101f13de31a3cec77653fca9a1 100644
--- a/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc
@@ -223,7 +223,7 @@ namespace Cyber{
   double Calo3DCluster::getEnergy() const{
     double result = 0;
     for(int m=0; m<m_2dclusters.size(); m++)
-      result += m_2dclusters[m]->getEnergy();
+      result += m_2dclusters[m]->getEnergy()*Escale;
     if(result == 0) result = getLongiE();
     return result;
   }
diff --git a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
index 8ae7e8a202e679a806243ef3763ecd60590a113b..dd53b6fe03b6f817ad9e42ff95d98958985b332f 100644
--- a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
+++ b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
@@ -11,12 +11,16 @@ namespace Cyber{
 
 
   StatusCode OutputCreator::CreateOutputCollections( CyberDataCol& m_DataCol,
+                                                     std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol,
                                                      DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler,
                                                      DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler,
                                                      DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler,
                                                      DataHandle<edm4hep::TrackCollection>& m_outTrkHandler,
                                                      std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler,
-                                                     DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler ){
+                                                     DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler,
+                                                     DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>& m_outMCPHitLinkHandler, 
+                                                     DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection>& m_outMCPClusterLinkHandler, 
+                                                     DataHandle<edm4hep::MCRecoParticleAssociationCollection>& m_outMCPPFOLinkHandler ){
 
     //Register the collections
     edm4hep::CalorimeterHitCollection* m_calohitCol = m_outRecHitsHandler.createAndPut();
@@ -27,11 +31,25 @@ namespace Cyber{
     edm4hep::ClusterCollection* m_Ecalcore    = m_outClusterColHandler["EcalCore"]->createAndPut();
     edm4hep::ClusterCollection* m_Hcalcluster = m_outClusterColHandler["HcalCluster"]->createAndPut();
 
+    //edm4hep::MCRecoCaloParticleAssociationCollection* m_mcphitLinkCol = m_outMCPHitLinkHandler.createAndPut();
+    edm4hep::MCRecoClusterParticleAssociationCollection* m_mcpclusterLinkCol = m_outMCPClusterLinkHandler.createAndPut();
+    edm4hep::MCRecoParticleAssociationCollection* m_mcppfoLinkCol = m_outMCPPFOLinkHandler.createAndPut();
+
+
     edm4hep::TrackCollection* m_trkCol = nullptr;
     if(settings.map_boolPars.at("UseTruthTrk") || settings.map_boolPars.at("UseTruthMatchTrk") ) m_trkCol = m_outTrkHandler.createAndPut();
 
 
-    //PFO
+    //Readin the MCP-CaloHit link into one collection
+    std::vector<edm4hep::MCRecoCaloParticleAssociation> m_readinmcphitLinkCol;
+    for(unsigned int ihit=0; ihit<settings.map_stringVecPars.at("CaloHitCollections").size(); ihit++){
+      if( map_CaloParticleAssoCol.find(settings.map_stringVecPars.at("CaloHitCollections")[ihit])!=map_CaloParticleAssoCol.end() ){
+        const edm4hep::MCRecoCaloParticleAssociationCollection* const_MCPCaloAssoCol = map_CaloParticleAssoCol[settings.map_stringVecPars.at("CaloHitCollections")[ihit]]->get();
+        for(int ilink=0; ilink<const_MCPCaloAssoCol->size(); ilink++) m_readinmcphitLinkCol.push_back(const_MCPCaloAssoCol->at(ilink));
+      }
+    }
+
+    //Create PFO
     std::vector<std::shared_ptr<Cyber::PFObject>> p_pfos = m_DataCol.map_PFObjects[settings.map_stringPars.at("OutputPFO")];
 
     for(int ip=0; ip<p_pfos.size(); ip++){
@@ -39,6 +57,7 @@ namespace Cyber{
       std::vector<const Track*> vec_trks = p_pfos[ip]->getTracks();
       std::vector<const Calo3DCluster*> vec_Ecalclus = p_pfos[ip]->getECALClusters();
       std::vector<const Calo3DCluster*> vec_Hcalclus = p_pfos[ip]->getHCALClusters();
+      std::map<const edm4hep::MCParticle, float> mc_PFOmap; mc_PFOmap.clear();
 
       double ecalcalib = 1.;
       double hcalcalib = vec_trks.size()==0 ? settings.map_floatPars.at("HCALNeutralCalib") : settings.map_floatPars.at("HCALChargedCalib");
@@ -68,9 +87,10 @@ namespace Cyber{
    
           auto _hit = m_calohitCol->create();
           _hit.setCellID(_cellID);
-          _hit.setEnergy( p_hit->getEnergy()*ecalcalib );
+          _hit.setEnergy( p_hit->getEnergy()*p_clus->getEnergyScale()*ecalcalib );
           _hit.setPosition( pos );
           _hit.setType(1); //Ecal barrel
+
           m_clus.addToHits(_hit);
         }
    
@@ -107,6 +127,19 @@ namespace Cyber{
         m_clus.setEnergy( tmp_clusE );
         m_clus.setPosition( pos );
         m_clus.addToClusters(m_core);
+
+        std::vector< std::pair<edm4hep::MCParticle, float> > mc_map = p_clus->getLinkedMCP();
+        for(auto iter : mc_map){
+          auto rel_MCP = m_mcpclusterLinkCol->create();
+          rel_MCP.setSim(iter.first);
+          rel_MCP.setRec(m_clus);
+          rel_MCP.setWeight( iter.second );
+
+          if( mc_PFOmap.find(iter.first)==mc_PFOmap.end() ) mc_PFOmap[iter.first] = m_clus.getEnergy()*iter.second;
+          else mc_PFOmap[iter.first] += m_clus.getEnergy()*iter.second;
+
+        }
+
         m_pfo.addToClusters( m_clus );
 
         EcalClusE += tmp_clusE;
@@ -119,7 +152,7 @@ namespace Cyber{
       for(int ic=0; ic<vec_Hcalclus.size(); ic++){
         auto p_clus = vec_Hcalclus[ic];
         auto m_clus = m_Hcalcluster->create();
-   
+
         for(int ih=0; ih<p_clus->getCaloHits().size(); ih++){
           const Cyber::CaloHit* p_hit = p_clus->getCaloHits()[ih];
           //auto _hit = m_hcalhitCol->create();
@@ -132,6 +165,17 @@ namespace Cyber{
         m_clus.setEnergy( tmp_clusE );
         edm4hep::Vector3f pos( p_clus->getHitCenter().x(), p_clus->getHitCenter().y(), p_clus->getHitCenter().z() );
         m_clus.setPosition( pos );
+
+        std::vector< std::pair<edm4hep::MCParticle, float> > mc_map = p_clus->getLinkedMCP();
+        for(auto iter : mc_map){
+          auto rel_MCP = m_mcpclusterLinkCol->create();
+          rel_MCP.setSim(iter.first);
+          rel_MCP.setRec(m_clus);
+          rel_MCP.setWeight( iter.second );
+
+          if( mc_PFOmap.find(iter.first)==mc_PFOmap.end() ) mc_PFOmap[iter.first] = m_clus.getEnergy()*iter.second;
+          else mc_PFOmap[iter.first] += m_clus.getEnergy()*iter.second;
+        }
         m_pfo.addToClusters( m_clus );
 
         HcalClusE += tmp_clusE;
@@ -144,6 +188,13 @@ namespace Cyber{
         TVector3 p3vec = vec_Pos*(  (EcalClusE+HcalClusE)/vec_Pos.Mag() );
         edm4hep::Vector3f p3(p3vec.x(), p3vec.y(), p3vec.z());
 
+        for(auto iter : mc_PFOmap){
+          auto rel_MCP = m_mcppfoLinkCol->create();
+          rel_MCP.setSim(iter.first);
+          rel_MCP.setRec(m_pfo);
+          rel_MCP.setWeight(iter.second/(EcalClusE+HcalClusE));
+        }
+
         m_pfo.setEnergy( EcalClusE+HcalClusE );
         m_pfo.setCharge( 0 );
         m_pfo.setMomentum( p3 );
@@ -170,24 +221,40 @@ namespace Cyber{
           m_pfo.setMomentum(p3);
           m_pfo.setMass(0.139570); //TODO: all charged particles are set to pion mass. 
           m_pfo.setEnergy( sqrt(vec_trks[trkIndex]->getMomentum()*vec_trks[trkIndex]->getMomentum() + 0.139570*0.139570) );
+
+          std::vector< std::pair<edm4hep::MCParticle, float> > mc_trkmap = vec_trks[trkIndex]->getLinkedMCP();
+
+          for(auto iter : mc_trkmap){
+            auto rel_MCP = m_mcppfoLinkCol->create();
+            rel_MCP.setSim(iter.first);
+            rel_MCP.setRec(m_pfo);
+            rel_MCP.setWeight(iter.second);
+          }
+
         }
         else{
           TVector3 p3vec = vec_Pos*(  (EcalClusE+HcalClusE)/vec_Pos.Mag() );
           edm4hep::Vector3f p3(p3vec.x(), p3vec.y(), p3vec.z());
 
+          for(auto iter : mc_PFOmap){
+            auto rel_MCP = m_mcppfoLinkCol->create();
+            rel_MCP.setSim(iter.first);
+            rel_MCP.setRec(m_pfo);
+            rel_MCP.setWeight(iter.second/(EcalClusE+HcalClusE));
+          }
+
           m_pfo.setEnergy( EcalClusE+HcalClusE );
           m_pfo.setCharge( 0 );
           m_pfo.setMomentum( p3 );
         }
 
       }
-//printf("  Create PFO #%d: charge %.1f, p4 (%.7f, %.7f, %.7f, %.7f), mass %.3f \n", ip, m_pfo.getCharge(), m_pfo.getMomentum().x, m_pfo.getMomentum().y, m_pfo.getMomentum().z, m_pfo.getEnergy(), m_pfo.getMass() );
+      //printf("  Create PFO #%d: charge %.1f, p4 (%.7f, %.7f, %.7f, %.7f), mass %.3f \n", ip, m_pfo.getCharge(), m_pfo.getMomentum().x, m_pfo.getMomentum().y, m_pfo.getMomentum().z, m_pfo.getEnergy(), m_pfo.getMass() );
     }
     return StatusCode::SUCCESS;
   }
 
 
-
   edm4hep::Track OutputCreator::TruthTrack(edm4hep::MCParticle _mcp, edm4hep::TrackCollection* _trkCol ){
     
     auto m_track = _trkCol->create();
diff --git a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
index f105c6c74aa0f88dc28a84345a293167fde29827..31b1c67642b201b93b81c6602ee72ff9751ec18b 100644
--- a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
+++ b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
@@ -65,7 +65,7 @@ namespace Cyber{
 
         for(int ilink=0; ilink<const_MCPTrkAssoCol->size(); ilink++){
           if( const_TrkCol[itrk] == const_MCPTrkAssoCol->at(ilink).getRec() ) {
-            m_trk->addLinkedMCP( std::make_pair(const_MCPTrkAssoCol->at(ilink).getSim(), const_MCPTrkAssoCol->at(ilink).getWeight()) );
+            m_trk->addLinkedMCP( std::make_pair(const_MCPTrkAssoCol->at(ilink).getSim(), const_MCPTrkAssoCol->at(ilink).getWeight()/(double)const_TrkCol[itrk].trackerHits_size() ) );
             break;
           }
         }